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NOMENCLATURE 


I . SYMBOLS 


Text 

a 


b 

BH 


B 


Z 


c 


Cf 


★ 

c 


If 


f 


w 


Computer 

Program 


FAH 

FAZ 

B 

FBH 

FBZ 

CF 


Definition 


Coefficient  in  the  trial  tempera- 
ture profile,  defined  by  Eq.  (19) 

Coefficient  in  solution  to  mechani- 
cal energy  equation,  defined  by 
Eq.  (16) 

Coefficient  in  solution  to  momen- 
tum equation,  defined  by  Eq . (15) 

Coefficient  in  the  trial  tempera- 
ture profile,  defined  by  Eq.  (20) 

Coefficient  in  solution  to  mechani- 
cal energy  equation,  defined  by 
Eq.  (16) 

Coefficient  in  solution  to  momen- 
tum equation,  defined  by  Eq . (15) 

Coefficient  in  the  trial  tempera- 
ture profile,  defined  by  Eq.  (21) 

Local  skin  friction  coefficient. 


f 2 

Sonic  speed  at  M = 1 

Second  derivative  of  dimensionless 
stream  function  used  by  Smith  and 
Clutter  [26 ] , 


F. ,F_ ,F_ , FF1,FF2,FF3,  Universal  functions  in  Eqs . (5), 
Fj,F^,Fg  FF4,FF5,FF6  (7),  and  (31) 


Blanks  indicate  that  additional  alphanumeric  symbols  may  be 
added  for  further  identification,  e.g.,  corresponding  to  sub- 
script notation. 


VI 1 


Text 

Computer 

Program 

Definition 

H 

H 

Incompressible,  adiabatic  shape 
factor, 

. - (Sj)" 

" WZ 

* 

H 

H S 

Shape  factor, 

* 63 

H = *1 

IDENT 

Computer  variable  used  at  Stan- 
ford conference  [18]  for  flow 
identification  purposes 

k 

Thermal  conductivity 

K 

CHT 

Heat  transfer  correction  para- 
meter introduced  in  Eqs . (30) 
and  (31) 

L 

Length  measured  along  body  sur- 
face, Eqs.  (37)  and  (38) 

[L] 

[L] 

Denotes  dimensions  of  length 

M 

M 

Mach  number 

* 

M 

M S 

Dimensionless  speed  ratio, 

MW 

Molecular  weight 

n 

N 

Exponent  of  Rg _ in  definition  of 

Z,  equal  to  1 for  laminar  boundary 
layers,  0.268  for  turbulent  boun- 
dary layers 

P 

Pressure 

Pr 

Prandtl  number 

q 

local  heat  flux 

Q 

Total  heat  transferred  from  wall 
to  boundary  layer,  defined  in  Eqs. 

(37)  and  (38) 

R 

R 

Body  radius 

Vlll 


Text 


Computer 

Program 


Definition 


r 

R«S  2 

Reoo 


]R 

s 

St 

T 

TI 

u 

w 

x,X 

Y/Y 

z 

a 


RL,  RT 

RD  2 


RE  INFL 


SL,ST 

T 

FST1NT 

W 

X 

Y 

Z 


Recovery  factor,  either  laminar 
or  turbulent 


Re>  idds  number  based  on  momentum 
thickness , 


Rfi  o 


P6u662 


w 


Unit  Reynolds  number  of  approach 
flow. 


Repo  __  PcqUqq 
L poo 

Universal  gas  constant 

Reynolds  analogy  factor,  either 
laminar  or  turbulent 


Stanton  number 


Temperature 

Freestream  turbulence  intensity, 
in  percent 

Streamwise  velocity  component 

Exponent  in  viscosity  power  law 
relation 

Streamwise  coordinate  defined  in 
Fig.  1;  in  Section  II  only,  x 
refers  to  coordinate  along  boun- 
dary 

Normal  coordinate  defined  in  Fig. 
1;  in  Section  II  only,  y refers 
to  coordinate  normal  to  boundary 

Momentum  parameter  defined  by 

Z = <S2R62n 

Angle  measured  from  stagnation 
point  on  circular  cylinder,  Fig. 
16 


IX 


Text 

3 


3 


E 


Computer 

Program  Definition 


BSTART  3tt  = BSTART *P I is  included  angle 

at  leading  edge  for  plane  or  axi- 
symmetric  external  flows.  Fig.  1 

Transformed  angle  for  axisymmetric 
flows , 

8 = £ 

PE  3-3 


Y 

6 


9 


IT 

A 


U 


G 

D 

D 1 


Specific  heat  ratio  of  ideal  gases 
Boundary  layer  thickness 
Displacement  thickness, 

6i = f(1  - f- 

1 o 1 p6  u6 


D 2 


Momentum  thickness, 


--D3 


Mechanical  energy  loss  thickness, 

63  = Jr*  r-f1  - tS->2idy 

J 0^6  6 u6 


04- 


Density  loss  thickness, 

6<  -Jh  ^ - 1,dy 


THET 


Temperature  ratio, 


6 = 


T - T 
e w 

T - T. 

e 6 


PI 

PPAR 


Constant,  3.14159 

Pohlhausen  parameter. 


A = 


Wu  du6 

yw  dx 


MU 


Absolute  viscosity 


x 


Computer 

Text  Program 


Definition 

Dummy  integration  variable  in 
Eqs . (37)  and  (38);  streamwise 

coordinate  along  boundary 

Density 

Shear  stress 


!• 


II.  SUBSCRIPTS  AND  SUPERSCRIPTS 


Text 

Computer 

Program 

Definition 

e 

Denotes  adiabatic  wall  conditions 

i , i+1, i-1 

1, 

IP1, 

IM1 

Subscripted  quantity  evaluated  at 
X = Xi'  Xi+1  °r  Xi-1 

inst 

INST 

Instability 

j f j+1 

Superscripted  quantity  evaluated 
at  (j)th  or  (j+1) st  iteration 

lam 

L 

Laminar 

m 

M 

Mean  value 

r 

Reference 

sep 

Separation 

trans 

TRANS 

Transition 

turb 

T 

Turbulent 

u 

— U 

Subscripted  quantity  depends  only 
on  velocity  profile 

w 

— w 

Subscripted  quantity  evaluated 
at  the  wall 

X 

Subscripted  quantity  based  on 
distance  from  stagnation  point 

6 

-D 

Subscripted  quantity  evaluated  at 
the  edge  of  the  boundary  layer 

o 

Refers  to  stagnation  conditions 

CO 

INF 

Refers  to  conditions  of  approach 
flow 

— 

B — 

Refers  to  average  values  over  in- 
tegration step 

Xll 


INTRODUCTION 


I . 


In  a large  number  of  flow  situations  of  engineering  interest, 
the  characteristics  of  the  developing  boundary  layer  are  of  prac- 
tical importance.  Examples  include  the  design  of  supersonic  noz- 
zles, flow  over  airfoils,  incipience  of  separation  in  pressure 
recovery  devices,  boundary  layer  growth  on  missile  bodies  and  af- 
terbodies and  the  resulting  effect  on  the  base  region,  etc. 

The  purpose  of  this  study  was  to  develop  an  easy-to-use 
FORTRAN  IV  computer  program  which  could  be  employed  routinely 
as  a boundary  layer  predictive  tool  for  a wide  variety  of  two- 
dimensional  ideal  gas  flows.  A number  of  the  available  methods 
were  considered,  with  the  rather  well-documented  [1-3]*  integral 
method  of  Walz  and  his  co-workers  being  chosen  for  use.  The  mo- 
tivation for  this  choice  was  in  part  due  to  the  fact  that  the 
method  described  in  [3]  was  one  of  the  four  integral  methods 
judged  "good"  at  the  Stanford  conference  on  incompressible  turbu- 
lent boundary  layers  and  in  part  due  to  the  flexibility  of  Compu- 
tational Method  II  which  is  detailed  in  [1].  This  method  may  be 
applied  to  incompressible  or  compressible,  laminar  or  turbulent, 
plane  two-dimensional  or  axisymmetric  boundary  layers  in  arbitrary 
pressure  gradients,  either  with  or  without  heat  transfer.  The 
computations  may  be  started  at  a stagnation  point,  a sharp  lead- 
ing edge,  or  at  an  arbitrary  streamwise  location  for  either  lami- 
nar or  turbulent  boundary  layers  with  the  proper  specification  of 
the  starting  parameters.  Transition  is  predicted  using  the  lami- 
nar instability  calculations  of  Wazzan,  et.  al.  [4],  together 
with  the  data  of  Granville  [5]  for  the  distance  between  the  in- 
stability and  transition  points.  Alternately,  the  transition 
point  can  be  specified  externally  as  an  input.  Separation  is 
also  predicted. 

Included  in  this  report  is  an  outline  of  the  theory  upon 
which  the  method  is  based,  a description  of  the  computer  program 
(COMPBL)  which  has  been  developed,  detailed  input  instructions 
and  example  input  and  output,  and  comparisons  of  computed  re- 
sults to  both  experiments  and  analytical  solutions  for  a number 
of  flow  fields.  A complete  program  listing  and  explanation  of 
error  messages  are  included  in  the  appendices. 


♦Numbers  in  brackets  refer  to  entries  in  REFERENCES. 
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II.  THEORETICAL  DEVELOPMENT 


A . General 

Computational  Method  II  of  Walz  [1],  on  which  program  COMPBL 
is  based,  utilizes  the  integral  momentum  and  mechanical  energy 
equations  for  the  boundary  layer: 


i du* 

— -3—  [2  + 

u dx  1 


-1  m 21 
62  ~ M<S 


d<53  , P 1 du6  r,  , 064  „ 2 

d^T  + 63'IT7  diT  [3  + 267  - M6  ] - 


2 

3 

itu,  i 

6 6 o 


T du  = 0 - (2) 


Momentum  equation  (1)  is  transformed  by  introducing  the  new 
dependent  variable, 


where 


to 


(4) 

(5) 


Fi  and  F2  are  universal  functions  which  may  be  evaluated  once  the 
trial  solutions  for  the  velocity  profiles  for  laminar  and  turbu- 
lent boundary  layers  and  the  coupling  law  between  the  tempera- 
ture and  velocity  profiles  are  specified. 


The  coordinate  x used  in  this  section  is  the  coordinate  along 
the  boundary  and  is  not  necessarily  the  same  as  the  coordinate 
x used  in  Section  III  or  program  COMPBL. 


3 


The  ratio,  ^ 

Pu  r i , u . 2,  , 

6.,  J P.u.  ^ u^  ^ 

H*  = / - — ^ ' (6) 

2 f Pu  ,,  u , j 

vr:  (1  - 5T>  ^ 

0 0 0 0 

is  introduced  into  the  mechanical  energy  equation  (2)  to  obtain 


H*  + H* 


F - — = 0 
3 Z 


where  F3  and  F4  are  also  universal  functions  dependent  on  the 
trial  velocity  profiles,  the  temperature-velocity  coupling  law, 
and  for  the  turbulent  case,  empirical  shear  stress  and  dissipa- 
tion integral  relations. 

Except  for  the  small  influence  of  R52  on  F4,  the  universal 
functions  may  be  reduced  to  depend  only  on  H,  M5 , and  0 where 
(x)  and  the  heat  transfer  parameter, 


0 ( x ) = 


T - T 
e w 

T - Tx 
e 6 


are  specified  inputs.  H is  defined  by: 

r<5 


H = (H* ) = 


( 6 , ) 

3 u 

(S0) 

2 u 


f u - <;r>2Hy 
U6  U6 


where  the  u subscript  denotes  that  the  enclosed  quantity  depends 
only  on  the  velocity  distribution  in  the  boundary  layer  (i.e., 
incompressible  and  adiabatic) . Noting  the  integral  definitions 
of  H*  and  H in  Eqs.  (6)  and  (9)  and,  in  particular,  the  similar 
manner  in  which  the  factor  P/P5  enters  the  integrands  in  the 
numerator  and  denominator  of  H*,  it  is  reasonable  to  expect  as 
a first  approximation  that. 


H*  = H. 


4 


The  final  relation  used  for 


H*  = H*(H,M6,6)  (11) 


is  based  on  the  consideration  of  limiting  cases  as  established 
by  Jischa  [6].  For  M5  = 0 this  relation  yields  H*/H  = 1 and  for 
Mfl  >>  1,  H*/H  - 1.18  (depending  on  the  value  of  0);  thus  the 
approximation  (10)  is  confirmed.  The  reduction  of  the  universal 
functions  for  the  general  compressible  case  with  heat  transfer 
to  H-dependent  functions  which  characterize  incompressible, 
adiabatic  (constant  property)  flows  requires  the  development  of 
additional  relationships  which  will  not  be  detailed  here.  The 
interested  reader  is  referred  to  [1] . 

Once  the  universal  functions  have  been  specified  as  func- 
tions of  H,  M5 , and  6 for  laminar  and  turbulent  boundary  layers, 
the  solution  proceeds  as  a step-by-step  simultaneous  integra- 
tion of  Eqs . (5)  and  (7),  with  Eq.  (11)  also  utilized.  This 

procedure  is  carried  out  with  a predictor-corrector  iterative 
scheme  as  follows.  Given  values  of  Hp  and  Zp  at  point  Xp  in 
the  boundary  layer,  a value  of  Hp+p  at  Xi+i  is  first  predicted 
(using  dH/dx  between  xp  and  Xi_i) . Average  values  of  the  uni- 
versal functions  over  the_integration  step  Ax  = (xp+p  - xp)  can 
then  be  evaluated,  e.g.,  Fp  = Fp((Hp  + Hp+p)/2,  (M£p  + M$p+p)/2, 
(0i  + 0p+p ) /2 ).  Using  these  average  values  and  approximating  the 
variation  of  the  freestream  velocity,  u<5,  as  a straight  line 
over  the  interval  Ax,  i.e.. 


u6(x)  = (u6).  + 
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Eqs.  (5)  and  (7)  become  linear,  first-order,  ordinary  differen- 
tial equations.  They  are  integrated,  respectively,  to: 


and 
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The  value  of  at  x^+j  is  then  found  by  inverting  Eq.  (11) 

using  the  value  of  H^+i*  from  (14) . If  this  corrected  value  of 
+ agrees  with  the  predicted  value  to  within  an  amount  speci- 
fied by  a convergence  criterion,  the  boundary  layer  parameters 
such  as  displacement  and  momentum  thicknesses,  skin  friction 
coefficient,  etc.  are  determined.  If  not,  iterations  proceed 
until  convergence  for  Hi+i  is  obtained.  Note  that  this  scheme 
obviates  numerical  differentiation  of  the  given  freestream  velo- 
city distribution  and  that  the  justification  for  linearizing 
U5 (x)  over  the  interval  Ax  rests  on  taking  a large  number  of 
base  points,  x^  , since  any  function  can  be  approximated  as 
piecewise  linear  if  the  intervals  are  sufficiently  small. 


B . Laminar  Boundary  Layers 

For  laminar  boundary  layers,  the  universal  functions  are 
evaluated  using  a one-parameter  trial  solution  for  the  velocity 
distribution  based  on  the  Hartree  wedge-flow  profiles.  The 
separation  profile  is  given  when  H falls  to  a value  of  1.515 
in  which  case  the  skin  friction  simultaneously  vanishes,  viz., 


(H  ) n = 1.515. 
sep  lam 


(17) 


In  accelerated  flows  values  of  H up  to  1.7  are  encountered. 

The  exponent  n in  the  definition  of  Z,  Eq.  (3),  has  a value 
of  1.0  for  laminar  flow. 

The  development  of  the  universal  functions  for  both  laminar 
and  turbulent  boundary  layers  is  also  dependent  on  the  coupling 
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and  r is  the  recovery  factor  in  Eqs.  (19-21). 

Equation  (18)  is  derived  strictly  under  the  assumptions  that  Pr 
- 1 and  that  the  streamwise  wall  temperature  and  pressure  gra- 
dients vanish  (i.e.,  dTw/dx  = dp/dx  = 0).  However,  as  pointed 
out  by  White  [8],  this  approximation  is  "surprisingly  accurate 
for  gases  under  general  boundary-layer  flow  conditions  at  high 
and  low  Mach  numbers  and  for  laminar  and  turbulent  flow." 

Based  on  the  temperature  distribution  of  Eq . (18),  the 

dimensionless  local  heat  flux  is  given  by  the  Reynolds  analog 
expression 


where  s is  the  Reynolds  analogy  factor.  Section  E will  give 
an  alternate  method  of  calculating  the  heat  transfer  when  strong 
wall  temperature  and/or  pressure  gradients  are  present. 
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C . Turbulent  Boundary  Layers 

The  trial  solution  for  the  velocity  profile  used  in  develop- 
ing the  universal  functions  for  the  turbulent  boundary  layer  is 
a combination  of  the  turbulent  log  law  and  the  law  of  the  wake 
suggested  by  Coles  [9] . In  addition,  the  empirical  shear  stress 
relation  of  Ludwieg  and  Tillmann  [10]  and  the  empirical  dissi- 
pation laws  of  both  Rotta  [11 ] -Truckenbrodt  [12]  and  Felsch 
[13,1,2,3]  are  used.  The  Rotta-Truckenbrodt  relation  is  em- 
ployed for  flows  with  favorable  pressure  gradients,  while  in 
regions  of  adverse  pressure  gradient,  the  Felsch  law  is  used. 

The  temperature  and  velocity  profiles  are  again  coupled  using 
Eq.  (18). 

Near  the  separation  point  for  turbulent  boundary  layers, 
it  has  been  noted  that  the  velocity  profiles  actually  depend 
on  two  parameters,  rather  than  the  single  one  used  in  this 
analysis.  Hence,  there  is  no  single  value  of  H which  can  be 
assigned  to  the  turbulent  separation  point,  i.e.,  1.50  < 

(Hsep) turb  < 1-57.  However,  the  following  approximation  is 
used : 


(H  ).  . a (H  ).  = 1.515.  (23) 

sep  turb  sep  lam 

The  upper  bound  on  H for  turbulent  boundary  layers  is  2.0.  The 
Ludwieg-Tillmann  exponent,  n = 0.268,  is  used  in  the  defini- 
tion of  Z , Eq.  (3) . 


D.  Transition 


Transition  from  a laminar  to  a turbulent  boundary  layer 
is  accomplished  in  program  COMPBL  either  by  specifying  the 
axial  location  of  transition  as  an  input  (variable  XTRAN)  or 
by  normal  calling  of  the  transition  subroutine  (TRANS).  This 
subroutine  utilizes  a correlation  to  the  incompressible,  adia- 
batic, laminar  stability  computations  of  Wazzan,  et.  al.  [4] 
as  listed  in  White  [8].  The  instability  point  is  located  when 
(R(52)  , which  is  monitored  at  each  axial  location  in  the  laminar 
boundary  layer,  exceeds  (R62)inst  from  the  Wazzan  correlation. 
The  transition  location  is  then  found  by  calculating  the  mean 
Pohlhausen  parameter,  Am,  from: 


^ rtrans 

^xtrans  - xinst^  j. 

xinst 


A (x) dx 


(24) 
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at  each  streamwise  location  and  using  the  curve  fit  listed 
in  White  [8]  to  Granville's  [5]  data  for  the  distance  between 
the  instability  and  transition  locations.  In  addition,  a cor- 
rection factor  is  employed  to  account  for  the  effects  of  free- 
stream  turbulence  intensity  on  the  transition  location.  The 
expression  used  is  the  correlation  given  by  Korst,  et.  al. 

[14]  to  Granville's  flat  plate  data.  The  combination  of  these 
two  curve  fits  gives  the  following  result: 


(R(S2)trans  (R<S2)inst  + 

(450  + 400e60Xm) (900  - 760TI  + 159TI2)/900  (26) 


where  TI  is  the  freestream  turbulence  intensity,  in  percent. 
For  freestream  turbulence  intensities  exceeding  approximately 
2.16%,  the  laminar  instability  and  transition  points  coincide. 

Once  the  transition  point  has  been  located,  the  calcula- 
tions are  switched  from  the  laminar  to  the  turbulent  regime 
with  the  aid  of  the  following  relations: 


H,  , = H, 

turb  lam 

(27) 

and 

(62}turb  = (<52)lam’ 

(28) 

The  second  relation  implies  that. 


Zturb 


*2R6 


0.268 


Zlam 
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-0.732 

2 


(29) 


Thus  in  this  technique,  the  momentum  thickness,  62,  is  continu- 
ous across  the  transition  point,  whereas  other  of  the  derived 
boundary  layer  parameters  such  as  displacement  and  boundary 
layer  thicknesses,  skin  friction  coefficient,  etc.,  are  not, 
due  to  their  totally  different  behavior  for  laminar  and  tur- 
bulent boundary  layers. 
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E.  Heat  Transfer  Correction 


The  coupling  law  Eq.  (18)  was  derived  under  the  assump- 
tions that  the  Prandtl  number  is  approximately  unity  and  that  the 
pressure  and  streamwise  wall  temperature  gradients  vanish. 

In  the  situation  where  only  the  properties  of  the  velocity 
boundary  layer  are  of  predominant  interest,  this  relation  is 
adequate  even  if  the  assumptions  are  badly  violated,  since  in 
this  case  only  the  average  properties  of  the  temperature  pro- 
file enter  the  analysis.  But  this  is  not  the  case  when  the 
heat  transfer  is  desired;  then  the  actual  shape  of  the  tem- 
perature distribution  is  of  utmost  importance  since  the  local 
heat  flux  is  proportional  to  the  derivative  of  this  distri- 
bution in  the  direction  normal  to  the  wall. 

To  account  for  the  effects  that  the  pressure  and  stream- 
wise  wall  temperature  gradients  (with  Pr  = 1)  might  have  on 
the  calculation  of  the  heat  transfer,  Walz  [1]  substituted 
the  following  trial  solution  for  the  temperature  profile, 


= a + [b  + K(x)  ] £-  + [c  - K (x)  ] (-^-)2,  (30) 

u6  u6 

into  the  thermal  energy  integral  equation  to  obtain: 


K'  + KF5  - Fg  = 0.  (31) 


Note  that  Eq.  (18)  is  a special  case  of  Eq . (30)  with  K = 0. 
Part  of  the  rationale  behind  this  choice  is  that  analytical 
solutions,  such  as  [15],  have  shown  that  for  nonisothermal 
walls  the  local  heat  flux  need  not  vanish  when  the  temperature 
difference  (Te  - Tw)  vanishes  (i.e.,  b = 0)  as  would  be  the  case 
when  using  Eq.  (18) . To  account  for  this  fact,  Walz  in- 
cluded the  correction  factor  K(x)  in  the  coefficient  of  u/u<5 
in  the  trial  temperature  profile. 

The  coefficients  F5  and  F6  in  Eq.  (31)  are  again  uni- 
versal functions  and  are  evaluated  in  the  same  manner  as  ex- 
plained above  for  the  universal  functions  of  the  momentum  and 
mechanical  energy  equations.  Equation  (31)  is  solved  in 
parallel  with  Eqs.  (5)  and  (7),  so  that  the  complete  scheme 
involves  simultaneous  integration  of  the  momentum,  mechanical 
energy,  and  thermal  energy  integral  equations  together  with 
relation  (II)  for  H*(H,M6,0).  The  dimensionless  heat  flux 
in  this  case  is  given  by  a "modified  Reynolds  analogy": 


qw 

3 

p <5U6 


(32) 


As  shown  later,  however,  this  heat  transfer  correction  proce- 
dure has  met  with  only  moderate  success. 

The  theory  discussed  above  in  Sections  A,  B,  C,  and  E is 
a synopsis  of  that  developed  by  Walz  in  [1].  For  further  de- 
tails, this  reference  should  be  consulted. 
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III.  COMPUTER  PROGRAM 


A . General  Description 

Program  COMPBL  has  the  following  characteristics  and  capa- 
bilities: 

1.  Computational  Method  II  of  Walz  [1]  is  utilized  to  compute 
the  integral  parameters  for  two-dimensional  laminar  or 
turbulent  boundary  layers. 

2.  Arbitrary  pressure  gradients  are  handled  with  the  boundary 
layer  either  compressible  or  incompressible,  axisymmetric 
or  plane  two-dimensional,  and  with  or  without  heat  trans- 
fer . 

3.  The  transition  location  may  either  be  specified  as  an  in- 
put or  predicted  by  the  Wazzan  [4 ] -Granville  [5]  method 
described  in  Section  II-D. 

4.  Separation  is  predicted  when  the  value  of  the  shape  para- 
meter H = (63) u/ ( 6 2 ) u falls  to  1.515  in  which  case  the 
skin  friction  simultaneously  vanishes. 

5.  The  computations  may  be  started  with  the  boundary  layer 
either  laminar  or  turbulent  at  a stagnation  point,  a 
sharp  leading  edge,  or  an  arbitrary  axial  location. 

6.  Freestream  velocity  input  data  may  be  specified  by  the 
Mach  number  M or  the  dimensionless  velocity  ratio  M*  = 
u/c*,  while  wall  temperature  data  may  be  input  as  either 
e = (Te  - Tw)/(Te  - T 5)  or  Tw/Too- 

7.  Options  are  available  whereby: 

a.  Calling  of  the  transition  subroutine  is  suppressed 
for  boundary  layers  thought  to  remain  laminar; 

b.  The  heat  transfer  calculation  is  corrected  using  the 
method  described  in  Section  II-E;  and 

c.  Intermediate  values  of  H and  other  variables  are 
printed  for  debugging  purposes. 

8.  COMPBL  is  written  in  FORTRAN  IV  with  typical  run  times  on 
the  CDC  CYBER  175  of  h - 1 CPU  seconds. 

Following  is  a list  of  the  limitations  of  the  method  and  program: 

1.  The  program  is  limited  to  ideal  gas  flows  with  certain  of 
the  input  variables  defaulted  to  values  for  air.  However, 
these  defaults  are  easily  overridden  on  input  (see  Section 
III-B) . 

2.  The  effects  of  surface  roughness,  freestream  turbulence 
level  (except  for  its  influence  on  the  transition  loca- 
tion) , shock-boundary  layer  interactions,  and  varying 
stagnation  pressure  are  not  considered. 

3.  Transition  is  assumed  to  occur  at  a point  in  an  essentially 
discontinuous  manner. 

4.  The  phenomena  of  relaminarization  and  laminar  separation 
with  reattachment  are  not  treated  automatically.  These 
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cases  may  be  calculated,  however,  by  restarting  the  com- 
putations at  the  proper  location. 

5.  The  output  consists  of  only  the  integral  boundary  layer 
parameters.  Velocity  and  temperature  profiles  are  not 
calculated . 


A list  of  the  input  variables  read  from  file  INPUT  follows. 
The  first  four  (4)  variables  are  literal  variables: 


FLOW literal  variable  describing  the  initial  flow  regime 

of  the  boundary  layer;  equal  to  "LAMINAR"  or 
"TURBULENT" 

GEOM literal  variable  describing  the  two-dimensional 


geometry;  equal  to  "PLANE  2-D"  or  "AXISYM" 

MT YPE ....  literal  variable  describing  whether  velocity  input 
data  is  in  terms  of  Mach  number,  M,  or  M*  = u/c*; 
equal  to  "MACH"  or  "MSTAR" 

TT YPE . . . . 1 iteral  variable  describing  whether  wall  tempera- 
ture input  data  is  in  terms  of  0 = (Te  - Tw)/(Te~T5) 
or  Tw/Too,  equal  to  "THETA"  or  "TWTINF" 

The  next  seventeen  (17)  variables  are  entered  via  NAMELIST  BL: 

ZSTART ...  starting  value  of  Z (default  = 0.0);  see  Section 
III-C 

HSTART ...  starting  value  of  H (default  = 1.572);  see  Section 
III-C 

BSTART. . . (BSTART  * PI)  is  the  included  angle  at  the  leading 

edge  for  both  plane  two-dimensional  and  axisymmetric 
external  flows  (default  = 0.0);  see  Fig.  1 and  Sec- 
tion III-C 

MINF Mach  number  (or  M*  for  MTYPE  = "MSTAR")  of  approach- 

ing flow;  see  Fig.  1 

RE  I NFL ...  unit  Reynolds  number  of  approach  flow,  PooU^/yoo,-  see 
Fig . 1 

XTRAN....X  location  of  specified  transition  point  (default 

= 0.0) 

FST I NT . . . freestream  turbulence  intensity,  in  percent  (de- 


fault = 0.0);  used  only  in  transition  subroutine 
to  locate  transition  point 

G ratio  of  specific  heats  (default  = 1.405) 

W exponent  on  viscosity  power  law,  i.e.. 


rf-  = (^-)W  (default  = 0.7) 
' o o 


RL laminar  recovery  factor  (default  = 0.85) 

RT turbulent  recovery  factor  (default  = 0.88) 

SL laminar  Reynolds  analogy  factor  (default  = 0.80) 
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ST. 

EPS 


turbulent  Reynolds  analogy  factor  (default  = 0.82) 
convergence  criterion  variable,  e.g.,  convergence 
on  H if: 


H.j+1 


< EPS 


(default  = 1.0  E-4) 


NOTRAN ...  logical  variable  which  if  .TRUE,  suppresses  calling 
of  the  transition  subroutine  for  boundary  layers 
thought  to  remain  laminar  (default  = .FALSE.). 

Note  that  NOTRAN  = .TRUE,  is  equivalent  to  setting 
XTRAN  to  a value  greater  than  the  largest  X input 
location. 

HTCORR ...  logical  variable  which  if  .TRUE,  invokes  the  use 

of  the  heat  transfer  correction  procedure  described 
in  Section  II-E  for  calculation  of  the  local  dimen- 
sionless heat  flux  (default  = .FALSE.) 

ERROR ....  logical  variable  which  if  .TRUE,  causes  intermediate 
H values  and  variables  associated  with  the  turbu- 
lent dissipation  integral  and  heat  transfer  correc- 
tion parameter  to  be  printed  for  debugging  purposes 
(default  = .FALSE.) 


The  next  five  (5)  input  variables  give  the  local  informa- 
tion for  each  point  at  which  the  boundary  layer  parameters  are 
to  be  calculated: 


X axial  location  of  the  boundary  point;  see  Fig.  1 

Y normal  location  of  the  boundary  point;  see  Fig.  1 

M OR  M”.. local  freestream  Mach  number  or  m*  (depending  on 
the  value  of  MTYPE) 

R cross-sectional  radius  for  axisymmetric  bodies  or 

normal  distance  from  centerline  to  boundary  for 
plane  two-dimensional  bodies;  see  Fig.  1 


THETA  OR 

TWT I NF . . . local  wall  temperature  data  (depending  on  the  value 
of  TTYPE);  for  adiabatic  flows  THETA  = 0 

The  output  variables  written  to  file  OUTPUT  are: 


X 

Z 

H 


axial  location  of  the  boundary  point 
local  value  of  z = 62  Rfio11 
local  value  of 


H = 


(Vu 

<Vu 


f 5 — [1  - (— ) 2 ] dy 
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D 1 D 2 


D 1 

D 2 


D999 


u(y  = D999).  = 0>999 

u6  P6U6 

RD2 momentum  thickness  Reynolds  number,  R^  = — 

r 'Jw 
w 

CF local  skin  friction  coefficient,  c,  = ~ 

QD1M dimensionless  local  wall  heat  flux,  ^w 


„ 3 

P6U6 

The  variables  written  to  file  TAPE4  are: 


X axial  location  of  the  boundary  point 

RCORR corrected  radius  or  normal  location,  RCORR  = R + D1 


File  T APE 4 contains  corrected  boundary  coordinate  infor- 
mation consisting  of  the  axial  coordinate,  X,  and  the  corrected 
axisymmetric  radius  or  plane  two-dimensional  centerline  dis- 
tance obtained  by  adding  the  displacement  thickness  to  R. 

This  information  can  be  used,  for  example,  in  the  design  of  su- 
personic nozzles  where  the  wall  coordinates  obtained  from  the 
inviscid  solution  are  relieved  by  an  amount  equal  to  the  dis- 
placement thickness  or  for  external  flows  where  the  corrected 
body  shape  is  obtained  by  adding  the  displacement  thickness 
to  the  original  body  shape. 

A sketch  of  the  coordinate  system  setup  is  shown  in  Fig. 

1.  The  orientation  of  the  X-Y  axes  is  completely  arbitrary 
as  shown  by  the  possible  choices  X-Y,  X 1 -Y  ’ , and  X"-Y"  in  each 
of  the  two  cases  sketched.  The  only  exceptions  to  this  state- 
ment are  in  the  case  where  the  transition  location  is  speci- 
fied as  an  input.  In  this  situation  it  is  assumed  that  axial 
coordinate  X increases  in  the  downstream  direction  (which  is 
the  case  for  all  choices  shown  in  Fig.  1),  and  that  specified 
location  XTRAN  is  not  equal  to  0.0  (which  is  equivalent  to 
"OFF").  For  all  choices  of  the  X-Y  axes,  R remains  the  body 


shape  factor, 
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99.9%  boundary  layer  thickness,  i.e., 
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radius  for  axisymmetric  geometries  or  the  distance  between  the 
centerline  and  the  boundary  for  plane  two-dimensional  geome- 
tries. In  the  axisymmetric  case  R must  be  entered  while  for 
plane  bodies  R is  entered  only  if  the  corrected  boundary  co- 
ordinates are  desired. 

The  boundary  is  approximated  by  a series  of  straight- 
line  segments;  the  distance  between  any  two  input  points  is 
given  by: 


Ax 


[(Xr+l 


- x.)2  + 


(Yi+l 


- yL) 


2,h 


(33) 


For  increased  accuracy,  therefore,  use  of  a large  number  of 
base  points  is  recommended  since  both  the  boundary  location 
and  the  freestream  velocity  distribution  are  approximated  by 
piecewise  linear  functions.  The  program  also  tends  to  run 
faster  when  using  a larger  number  of  input  locations  because 
savings  in  the  number  of  iterations  required  for  convergence 
usually  outweighs  the  fact  that  more  points  are  calculated. 

It  should  also  be  noted  that  the  program  does  a certain  amount 
of  interval  subdivision  automatically  so  that  the  number  of 
base  points  actually  calculated  is  greater  than  the  number 
entered . 


Also  shown  in  Fig.  1 is  the  location  of  the  upstream 
variables  Moo(MINF),  Reoo/L(RE  INFL)  , and  Too . This  location  may 
be  either  the  undisturbed  approach  flow  as  shown  in  the  ex- 
ternal flow  of  Fig.  1(a)  or  the  freestream  adjacent  to  the 
first  input  point  as  shown  in  Fig.  1(b).  Too  is  the  approach 
flow  static  temperature  used  in  forming  the  dimensionless  wall 
temperature  ratio  Tw/T^  = TWTINF.  For  flows  in  which  the 
freestream  stagnation  conditions  are  known,  the  unit  approach 
Reynolds  number  Reoo/L  = RE  I NFL,  may  be  calculated  with  the 
aid  of  the  following  equation: 


RE INFL 


p M 
o 


00 


r — — 2.155  [i 
lMW  y j l 


Y + l 
2 (y-1) 


(34) 


The  meaning  of  variable  BSTART  is  also  shown  in  Fig.  1(a). 

The  input  variables  G,  W,  RL,  RT,  SL,  and  ST  have  been 
defaulted  to  values  corresponding  to  air.  For  gases  other  than 
air  the  appropriate  values  for  the  specific  heat  ratio,  G,  and 
the  viscosity  power  law  exponent,  W,  should  be  entered  with  the 
following  expressions  for  recovery  factors  and  Reynolds  analogy 
factors  being  recommended: 
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RL  = Pr1/2  RT  = Pr1/3  S L - ST  = Pr2,/3.  (35) 


For  each  base  point,  the  output  is  written  to  files  OUTPUT 
and  TAPE 4 and  is  identified  only  by  the  axial  coordinate,  X, 
of  the  point.  In  addition,  the  boundary  layer  parameters  are 
written  only  for  locations  corresponding  to  the  input  base 
points,  with  two  possible  exceptions.  If  ERROR  = .TRUE,  is 
specified,  the  output  at  all  of  the  calculated  points,  in- 
cluding the  automatic  subdivisions,  is  printed  together  with 
some  intermediate  results.  The  printout  may  become  quite 
long  in  this  case.  The  second  exception  is  that  the  output 
data  for  the  base  point  directly  downstream  from  a transition 
location  may  not  be  printed,  with  the  results  at  an  interme- 
diate subdivided  location  being  substituted.  This  is  a con- 
sequence of  the  automatic  interval  dividing  scheme  which  is 
used.  When  transition  is  to  be  predicted,  it  is  recommended 
that  an  initial  run  be  made  to  approximately  locate  the  laminar 
instability  and  transition  points  and  then  the  boundary  layer 
recalculated  with  a relatively  fine  gridpoint  spacing  around 
these  locations  so  that  the  instability  and  transition  points 
may  be  found  more  precisely. 

The  variables  Z and  H have  been  included  in  the  output 
to  provide  a means  for  monitoring  the  calculations  and  the 
state  of  the  boundary  layer.  As  mentioned  previously,  in  any 
region  where  H is  only  slightly  greater  than  1.515  for  laminar 
boundary  layers  or  where  1.515  < H < 1.57  for  turbulent  bound- 
ary layers,  the  onset  of  separation  is  imminent.  The  output 
values  of  Z and  H are  also  helpful  in  setting  the  values  of 
ZSTART  and  HSTART  when  the  boundary  layer  calculations  are 
restarted  at  some  downstream  location.  Downstream  restarting 
of  the  calculations  may  be  used,  for  example,  if  convergence 
problems  (hopefully  nonexistent)  are  encountered  in  the  in- 
tegration scheme  or  if  the  phenomena  of  relaminarization  or 
laminar  separation  with  reattachment,  which  are  not  handled 
automatically,  are  to  be  analyzed.  This  procedure  avoids  re- 
calculation of  the  upstream  portion  of  the  flow. 

The  local  dimensionless  heat  flux  which  is  written  to 
OUTPUT  is  defined  as: 


QDIM 


w 


P6U6 


(36) 


where  qw  is  the  local  wall  heat  flux.  The  total  heat  per  unit 
depth  transferred  from  the  wall  to  the  boundary  layer  for  plane 
two-dimensional  geometries  is  given  by: 


18 


(37) 


Q =J  qw(?)d?' 

while  for  axisymmetric  bodies  the  heat  transfer  is  given  by, 


Q = 2tt  J qw(5)RU)  d£  (38) 

where  £ is  measured  along  the  boundary.  Positive  signs  for 
qw  and  Q indicate  heat  transfer  from  the  wall  to  the  boundary 
layer. 

The  dimensional  variables  in  the  input/output  lists  are: 
ZSTART,  RE  I NFL,  XTRAN,  X,  Y,  R,  Z,  Dl,  D2,  D999,  and  RCORR . 

All  have  dimensions  of  length,  [L] , except  for  RE  I NFL  which 
has  dimensions  of  reciprocal  length,  [L]”f  The  only  require- 
ment that  need  be  met  concerning  units  is  that  input  variables 
ZSTART,  RE  I NFL,  XTRAN,  X,  Y,  and  R be  entered  with  the  same 
length  unit.  On  output  Z,  Dl,  D2  and  D999  will  have  this  same 
unit  of  length. 

Specific  instructions  regarding  the  input  procedure  to- 
gether with  an  example  are  given  in  Section  III-D. 


Starting  of  the  Computations 


In  order  to  begin  integrating  the  first  order  ordinary 
differential  equations  given  by  (5)  and  (7)  initial  values  of 
Z and  H must  be  specified.  This  is  accomplished  through  input 
variables  ZSTART  and  HSTART  as  follows. 


ZSTART  When  the  computations  are  started  at  a sharp  lead- 
ing edge  or  a stagnation  point,  the  value  ZSTART  = 0.0  is  used. 
If,  however,  the  calculations  are  started  at  an  arbitrary  stream- 
wise  location,  where  the  boundary  layer  thickness  is  not  zero, 
the  definition  of  Z is  employed: 

ZSTART  = _R<$  _n  n = 1.0  for  laminar  boundary 

layers 

(39) 

= 0.268  for  turbulent 
boundary  layers 


Hence,  the  momentum  thickness,  62,  must  be  known  at  the  start- 
ing location.  The  momentum  thickness  Reynolds  number,  R^/ 
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may  be  calculated  from  the  unit  Reynolds  number  of  the  approach 
flow,  Re  /L,  as: 


R, 


P6u6<52 
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(40) 


where  the  subscript  "i"  denotes  the  starting  location. 

HSTART  For  stagnation  points  or  sharp  leading  edges  the 
procedure  for  determining  HSTART  is  the  same  whether  the  bound- 
ary layer  is  started  in  the  laminar  or  turbulent  regimes  since 
the  first  step  is  always  calculated  as  laminar  (although  it 
can  be  made  arbitrarily  small) . For  plane  two-dimensional 
geometries  and  non-zero  included  wedge  angles,  Bn  = BSTART  :: 

PI  ^ 0,  the  starting  value  of  H in  both  the  incompressible 
and  compressible  cases  is  obtained  from  Table  I,  which  is  a 
compilation  of  the  Hartree  solutions,  since  in  the  immediate 
vicinity  of  the  leading  edge  the  flow  may  be  considered  incom- 
pressible. For  plane  flat  plate  flows,  Bit  = 0,  the  starting 
value  of  H is  obtained  from  Fig.  2,  which  details  the  effect 
of  compressibility  (Moo)  on  HSTART. 

For  axisymmetric  geometries  with  subsonic  approach  flow, 

0 < 1,  HSTART  is  obtained  from  Table  I using  the  trans- 

formed angle  Be  where: 


B 

E 3 - B ' 


(41) 


For  supersonic  flow,  Moo  > 1 , approaching  on  axisymmetric  body, 
the  flow  fields  along  the  cone  approximating  the  bow  of  the 
body  and  along  a flat  plate  are  similar  so  that  Fig.  2 is  also 
used  in  this  case. 

To  determine  HSTART  for  boundary  layer  calculations  ini- 
tiated at  an  arbitrary  axial  location,  either  the  local  skin 
friction  coefficient,  cf,  or  the  displacement  thickness,  6^, 
must  be  known  at  the  starting  location  in  addition  to  the 
momentum  thickness  which  is  required  for  ZSTART.  For  laminar 
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boundary  layers  HSTART  may  be  determined  by  solving  either 
of  the  following  transcendental  equations  for  H: 


Cf  known: 
or 

e 1 known : 

<50 


(H  - 1.515) 


0.7158 


cf  R«2 


[1+r^  M6i2  (H  - 0i) (2  - H) ] 3-4522 


(42) 


4.0306  - 4 . 2845  (H  - 1 . 515)  0 * 3886 ] [1  + r M.5..2 (H  - 6^(2  - H) ] 


+ r - e,)  - 


(43) 


For  turbulent  boundary  layers,  the  following  equations  are  used: 

0.268 

(44) 


known: 


(H  - 1.515) 
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R6- 


[1  + r^  M6;L2(H  - Of)  (2  - H)  ] 


0.07788 


or 


-r—  known: 
62 


[1  + 1.48  ( 2-H)  + 104  (2-H)6*7]  [1  + r ^ M^2  (H-Gf)  ( 2-H)  ] 

•.  2 <5  •> 

+ r M^tH-Gf)  = (45) 


The  approximation  H*  - H is  used  in  the  above  equations  and 
is  completely  adequate  for  the  purpose  of  determining  HSTART. 
For  incompressible  boundary  layers,  <<  1,  the  above  equa- 

tions are  greatly  simplified  and  all  but  (45)  can  be  solved 
explicitly  for  H.  For  cases  with  heat  transfer  where  the  Stan- 
ton number,  St,  is  known,  HSTART  may  be  found  from  Eqs.  (42) 
or  (44)  using  cf  = (St) (2s) , assuming  that  the  Reynolds  ana- 
logy is  accepted. 

BS TART  Variable  BSTART  is  required  only  for  boundary 
layer  calculations  started  at  sharp  leading  edges  or  stagna- 
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tion  points.  For  either  ax i symmetric  or  plane  geometries 
BSTART  ::  PI  = fin  is  the  total  included  angle  at  the  nose. 

A summary  of  the  procedure  used  in  determining  the  start- 
ing values  ZSTART,  HSTART,  and  BSTART  is  presented  in  Table 
II. 


D . Input  Instructions  and  Example 

The  input  variables  listed  in  Section  III-C  are  entered 
in  the  following  way.  The  first  card  (record)  is  used  to  en- 
ter a title  which  is  printed  on  files  OUTPUT  and  TAPE4  to  iden- 
tify the  results.  Any  message  up  to  80  columns  can  be  used, 
but  for  aesthetic  reasons  the  title  should  be  centered  in  the 
80  columns.  It  is  suggested  that  the  length  unit  for  dimen- 
sion [L]  be  entered  as  part  of  this  title.  On  the  next  card 
the  four  literal  variables  FLOW,  GEOM  MTYPE,  and  TTYPE  are 
input  in  4A10  format.  The  value  of  these  literals  should  be 
left  justified  in  each  10  column  block.  NAMELIST  "BL"  which 
encompasses  input  variables  ZSTART  - ERROR  is  input  via  the 
next  card(s).  The  first  column  on  each  NAMELIST  card  should 
be  blank.  All  variables  in  "BL"  except  MINF  and  REINFL  are  de- 
faulted so  that  at  a minimum  only  these  two  and  the  selected 
variables  to  be  overridden  need  be  entered.  The  remaining 
cards  (records)  contain  the  local  data:  X,  Y,  M (or  M::)  , R, 
THETA  (or  TWTINF)  for  each  point  at  which  the  boundary  layer 
parameters  are  to  be  found.  These  variables  are  entered  in 
5F10  format  and  as  many  points  as  desired  can  be  used.  For 
boundary  layers  starting  at  sharp  leading  edges  or  stagnation 
points,  the  first  of  these  local  data  cards  should  contain  the 
information  for  the  tip.  In  the  case  where  the  computations 
are  started  at  an  arbitrary  axial  location,  the  first  local 
data  card  should  consist  of  the  information  for  the  point  at 
which  ZSTART  and  HSTART  were  determined. 

As  an  example,  the  flow  over  a NACA  0012  airfoil  at  zero 
angle  of  attack  is  considered.  This  example  is  taken  from 
the  report  by  McNally  [16]  who  cites  Becker  [17]  as  the  ori- 
ginal source  of  the  data.  The  airfoil  profile  and  freestream 
M“  distribution  are  shown  in  Fig.  4(b). 

A schematic  of  the  input  file  for  this  case  is  shown  in 
Fig.  3.  The  first  card  contains  the  centered  title,  while 
the  second  contains  the  input  values  of  literal  variables 
FLOW,  GEOM,  MTYPE,  and  TTYPE.  For  this  example,  the  boundary 
layer  is  assumed  to  begin  in  the  LAMINAR  regime;  the  airfoil 
has  PLANE  2-D  geometry;  the  velocity  input  data  is  in  terms 
of  MSTAR;  and  the  wall  temperature  data  is  in  terms  of  TWTINF. 
Note  that  the  value  of  each  literal  is  left  justified  in  its 
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10  column  block.  NAMELIST  BL  requires  only  a single  card  for 
this  example  and  as  required  it  begins  with  the  first  column 
blank  and  the  characters  $BL  and  ends  with  $.  Since  the  bound- 
ary layer  calculations  are  started  at  a stagnation  point, 

ZSTART  = 0.,  but  this  is  the  default  value  so  it  isn't  entered. 
From  the  average  angle  between  the  first  two  input  points,  it 
is  easily  determined  that  BSTART  = 0.73951.  Entering  Table  I 
with  this  value  of  3 = BSTART  yields  HSTART  = 1.6190.  From 
the  test  conditions  described  in  [16],  the  following  values 
for  the  approach  flow  M::  and  unit  Reynolds  number  are  obtained: 
MINF  = 0.30863  and  RE  I NFL  = 2.04146  X 106  [ft“lj.  Since  a 
normal  call  to  the  transition  subroutine  is  to  be  made,  the 
default  values  of  both  XTRAN  and  NOTRAN  are  used.  Also  the 
freestream  turbulence  intensity  is  assumed  to  be  zero,  the 
flowing  gas  is  air,  and  neither  the  heat  transfer  correction 
nor  the  error  options  are  desired;  consequently,  all  other 
variables  except  G (y)  are  left  at  their  default  values.  G 
is  overridden  from  1.405  to  1.4.  The  remaining  cards  con- 
tain the  local  data  X,  Y,  MSTAR,  R,  and  TWTINF  for  each  base 
point  in  5F10  format.  Note  that  R values  are  not  necessary 
and  therefore  not  entered  for  this  plane  two-dimensional  geo- 
metry since  corrected  boundary  coordinate  data  is  not  desired 
in  this  case.  The  airfoil  is  assumed  to  be  isothermal  at  the 
approach  flow  stagnation  temperature,  T0  ^ so  that  Tw/T  OO  — 
TWTINF  = 1.0161.  In  this  example  the  unit  for  input  variables 
(REINFL)-l,  X,  and  Y is  feet. 

The  listing  of  file  OUTPUT  for  this  example  flow  is  shown 
in  Fig.  4(a).  The  locations  of  the  laminar  instability  and 
transition  points  are  clearly  indicated,  as  shown.  Fig.  4(h) 
compares  the  results  computed  from  COMPBL  with  the  experimental 
data . 

As  an  example  for  calculating  corrected  wall  coordinates, 
turbulent  boundary  layer  computations  were  carried  out  from 
the  throat  to  the  exit  of  a plane  two-dimensional  Mach  2 noz- 
zle. The  listing  of  file  TAPE4  for  this  example  is  presented 
in  Fig.  5. 
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IV.  COMPARISONS  TO  EXPERIMENT  AND  ANALYTICAL  SOLUTIONS 


In  order  to  test  the  accuracy  of  both  the  method  and  the 
computer  program,  boundary  layer  computations  were  carried 
out  for  a number  of  flow  situations  for  which  either  data  or 
analytical  solutions  are  available.  The  resulting  comparisons 
are  shown  in  Figs.  6-21.  In  addition,  the  results  of  calcu- 
lations for  a conical  boattail,  a cylindrical  missile  body/ 
conical  boattail  combination,  and  a plane  two-dimensional 
Mach  2 nozzle  are  presented  in  Figs.  22-24  to  demonstrate 
how  the  boundary  layer  computations  can  be  coupled  to  inviscid 
solutions.  Unless  otherwise  indicated,  adiabatic  wall  condi- 
tions (THETA  H 0.0)  are  used  and  the  ideal  gas  is  air. 


A . Turbulent  I n compres s i b 1 e Boundary  Layers 

Since  the  turbulent  regime  is  probably  more  common  and 
therefore  of  more  practical  importance  than  the  laminar,  a 
arge  number  of  calculations  were  made  for  this  case.  The 
five  example  flows  analyzed  for  the  incompressible  turbulent 
regime  were  all  test  cases  used  at  the  Stanford  conference 
and  were  taken  from  reference  [ 1 8 J . In  each  case  the  value 
of  HSTART  was  determined  by  taking  the  average  of  the  H values 
found  from  Eqs.  (44)  and  (45).  The  comparisons  of  the  compu- 
tations to  the  data  are  shown  in  Figs.  6-10. 

The  first  case  is  the  Ludwieg  and  Tillmann  accelerating 
flow  (I  DENT  = 1300) + which  may  be  considered  as  relatively 
"easy"  to  predict.  As  shown  in  Fig.  6,  the  agreement  between 
the  calculations  and  the  data  is  excellent.  The  remaining 
four  example  flows  can  be  placed  in  the  "difficult"  category, 
being  either  strong  adverse  pressure  gradient,  separating, 
or  relaxing  flows.  The  results  for  the  Ludwieg  and  Tillmann 
strong  adverse  pressure  gradient  case  (1200)  are  shown  in 
Fig.  7.  The  agreement  is  good  up  to  about  X = 3m  at  which 
point  the  data  tends  toward  separation  while  the  computations 
do  not.  In  this  region,  however.  Coles  and  Hirst  [18]  point 
out  that  the  data  contains  large  discrepancies  in  the  momen- 
tum balance.  Figure  8 presents  the  comparison  for  Clauser, 
flow  number  1 (2200) . The  computations  and  data  agree  reason- 
ably well  except  perhaps  in  the  downstream  section  where  62 
is  overestimated  and  cf  is  underestimated.  Large  momentum 
balance  discrepancies  are  again  noted  for  this  region.  The 
comparison  for  Moses,  case  3 (3800)  which  is  a case  of  ex- 
tremely strong  adverse  pressure  gradient  leading  to  separa- 
tion, is  shown  in  Fig.  9.  The  data  indicates  separation  at 


+Computer  variable  I DENT  was  used  at  the  Stanford  conference 
on  turbulent  boundary  layers  for  flow  identification  purposes. 
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approximately  X = 19.2  in.  while  the  calculations  do  not  in- 
dicate separation  at  all.  Although  momentum  balance  discre- 
pancies again  occur  for  the  data  in  the  downstream  section, 
they  probably  cannot  fully  explain  the  lack  of  agreement  be- 
tween the  data  and  predictions.  The  final  incompressible 
turbulent  boundary  layer  calculation  is  the  Bradshaw  and 
Ferriss  relaxing  flow  (2400)  shown  in  Fig.  10.  The  predic- 
tions and  data  are  in  good  agreement  except  for  the  underpre- 
diction of  the  skin  friction  coefficient,  cf,  at  the  down- 
stream stations. 

Of  all  the  flow  cases  used  for  the  Stanford  conference, 
Coles  and  Hirst  [18]  list  those  data  sets  which  they  feel  are 
preferable  because  of  experimental  reliability,  superior  in- 
strumentation, etc.  The  Ludwieg-Tillmann  accelerating  flow 
and  the  Bradshaw-Ferriss  relaxing  flow  fall  into  this  cate- 
gory, and  for  these  cases  the  program  results  agree  well  with 
the  data.  In  any  event,  since  the  nominal  flow  situation  is 
much  less  difficult  than  the  last  four  examples  considered 
above,  it  is  valid  to  conclude  that  the  program  developed 
herein  is  reasonably  accurate  for  incompressible  turbulent 
boundary  layers. 


B . Turbulent  Compressible  Boundary  Layers 

Five  test  boundary  layers  were  also  computed  for  the  com- 
pressible turbulent  regime.  The  first  of  these,  shown  in  Fig. 
11,  is  a very  strong  adverse  pressure  gradient  supersonic  flow 
reported  by  McLafferty  and  Barber  [19].  The  agreement  between 
the  calculations  and  the  data  is  relatively  good,  although 
the  streamwise  gradients  of  the  prediction  curves  are  not  as 
large  as  those  indicated  by  the  data.  Notice  that  the  bound- 
ary layer  is  very  near  separation  for  X > 2.5  in.  It  should 
be  mentioned  that  rather  than  determining  HSTART  from  Eq.  (45) 
for  the  given  initial  value  of  61/62,  HSTART  was  obtained  by 
computing  the  supersonic  flat  plate  boundary  layer  at  the  given 
entrance  conditions  in  the  streamwise  direction  until  the  ini- 
tial entrance  value  of  62  was  matched.  This  method  was  used 
since  it  modelled  well  the  entrance  region  to  the  test  section 
and  since  the  data  points  for  the  shape  factor,  61/62 , were 
taken  from  a small  graph  in  [19]. 

The  next  three  cases  are  taken  from  the  report  by  Winter, 
Smith,  and  Rotta  [20]  and  consist  of  the  flow  over  a waisted 
body  of  revolution  at  approach  Mach  numbers  of  Moo  = 0.6,  1.4, 
and  2.8.  The  measured  freestream  Mach  number  distributions 
indicate  that  the  pressure  gradient  is  initially  adverse  (for 
x/H  > 0.4)  followed  by  a relaxation  region.  The  average  value 
of  H determined  from  Eqs . (44)  and  (45)  was  used  for  HSTART. 
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As  shown  in  Figs.  12,  13,  and  14  the  shape  factor,  61/62,  is 
well  predicted  in  all  three  cases  as  is  the  momentum  thickness 
for  Moo  = 0.6  and  1.4.  However,  the  skin  friction  coefficient 
and  the  momentum  thickness  for  Moo  = 2.8  are  overpredicted, 
particularly  in  the  region  of  the  body  waist  (x/SL  = 0.7). 
Lewis,  Kubota,  and  Webb  (21)  have  noted  that  the  integral  mo- 
mentum balance  is  not  satisfied  in  this  region,  although  it 
is  doubtful  that  this  can  account  entirely  for  the  discrep- 
ancy between  the  data  and  the  computations. 

The  final  compressible  turbulent  flow  considered  is  the 
transonic  shock  wave/boundary  layer  interaction  case  investi- 
gated by  Alstatt  [22].  The  freestream  Mach  number  variation, 
shown  in  Fig.  15,  contains  regions  of  both  strong  accelera- 
tion and  strong  deceleration  leading  to  separation.  Since 
only  displacement  thickness  measurements  were  reported,  the 
values  of  ZSTART  and  HSTART  were  determined  by  computing  the 
flat  plate  boundary  layer  at  the  entrance  conditions  in  the 
streamwise  direction  until  the  entrance  value  of  6g,  was 
matched.  As  shown  in  Fig.  15,  the  measured  and  computed 
values  of  the  displacement  thickness  agree  very  well,  except 
perhaps  near  X = 14  in.  where  the  calculations  show  a sharp 
peak.  The  separation  point  was  computed  as  X = 21.4  in. 
agreeing  well  with  the  measured  value  of  X = 21.75  in. . 

Alstatt  [22]  reported  an  instability  problem  in  the  boundary 
layer  solution  technique  of  Nash  and  Hicks  [23]  for  a number 
of  different  eddy-viscosity  models  occurring  just  downstream 
of  X = 16  in.;  no  such  problems  were  encountered  here. 

Since  all  of  the  cases  discussed  above  for  the  compres- 
sible turbulent  regime  are  considered  as  "difficult"  and 
since  the  computations  and  data  agree  reasonably  well  for  the 
most  part,  it  can  be  concluded  that  program  COMPBL  yields  ac- 
curate approximations  to  the  boundary  layer  behavior  for  this 
regime . 


C . Laminar  Boundary  Layers 

For  laminar  boundary  layers,  three  examples  were  calcu- 
lated as  shown  in  Figs.  16,  17,  and  18.  All  three  boundary 
layers  begin  at  a plane  stagnation  point,  HSTART  = 1.625, 
BSTART  = 1.0.  The  first  two  cases  are  incompressible  non- 
similar flows  over  a circular  cylinder  investigated  by  Hie- 
menz  [24]  and  over  an  elliptic  cylinder  reported  by  Schubauer 
[25].  For  each  of  these  cases,  Smith  and  Clutter  [26]  have 
computed  the  corresponding  laminar  boundary  layers  using  their 
well  known  finite  difference  technique,  while  for  the  ellip- 
tic cylinder  Hartree  [27]  reported  earlier  hand-performed  cal- 
culations. In  both  cases,  the  agreement  between  the  results 
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computed  here  from  the  Walz  approximation  theory  and  the  "exact" 
finite  difference  calculations  is  remarkable.  The  comguted 
separation  point  for  the  circular  cylinder  is  a - 80.3°  agree- 
ing well  with  Hiemenz ' observed  value  of  a = 80°.  Neither 
the  present  results  nor  Smith  and  Clutter's  calculations  indi- 
cate separation  for  the  elliptic  cylinder,  although  Schubauer 
measured  it  to  occur  at  X = 1.99  +0.02.  However,  using  a 
polynomial  fit  to  Schubauer' s measured  pressure  distribution, 
rather  than  the  one  given  by  Hartree,  yielded  separation  at 
X = 1.96  with  the  present  method.  The  results  shown  in  Fig. 

17  are  for  the  Hartree  fit  to  the  pressure  distribution  since 
this  is  the  one  upon  which  the  Smith  and  Clutter  calculations 
were  based. 

Figure  18  compares  the  results  of  the  present  computa- 
tion scheme  with  the  finite  difference  calculations  of  Flugge- 
Lotz  and  Eichelbrenner  [28]  (taken  from  [29])  for  a laminar, 
compressible  airfoil-type  flow.  The  agreement  is  good  except 
that  the  separation  point  is  predicted  to  be  somewhat  too  far 
downstream. 

Based  on  these  examples,  the  present  method  appears  to 
be  very  accurate  for  laminar  boundary  layers. 


D . Boundary  Layers  with  Heat  Transfer 

Two  heat  transfer  cases  were  considered.  The  first,  shown 
in  Fig.  19,  consists  of  a laminar,  Mach  3 flow  over  a flat 
plate  with  the  wall  temperature  distribution  shown  in  the 
figure.  An  analytical  solution  of  this  problem  is  given  by 
Chapman  and  Rubesin  [15].  The  computation  of  the  heat  trans- 
fer based  on  the  simple  Reynolds  analogy,  K = 0,  is  greatly 
in  error  as  may  be  expected  for  the  strong  axial  wall  tempera- 
ture gradient  present  here.  Using  the  heat  transfer  correc- 
tion procedure  ("modified  Reynolds  analogy")  discussed  in  Sec- 
tion I I-E , K = K (x) , puts  the  computed  heat  transfer  into  much 
better  agreement  with  the  analytical  solution  although  it 
seemingly  overcorrects  for  the  effect  of  variable  wall  tem- 
perature . 

The  second  heat  transfer  case,  taken  from  the  report  by 
Boldman,  Schmidt,  and  Ehlers  [30],  is  a turbulent  boundary 
layer  flow  on  the  walls  of  a cooled,  supersonic  nozzle.  Both 
the  wall  temperature  and  freestream  Mach  number  vary  strongly 
in  the  axial  direction  as  shown  in  Fig.  20.  The  heat  transfer 
computations  based  on  the  Reynolds  analogy,  K = 0,  do  not  a- 
gree  well  with  the  data  and  err,  as  expected,  on  the  high  side 
for  the  strong  favorable  pressure  gradient  within  the  nozzle. 
However,  the  "corrected"  results,  which  are  not  presented,  are 
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in  even  poorer  agreement  with  the  data.  As  for  the  case  above, 
the  corrections  are  too  strong  and  in  fact  predict  a change  in 
the  direction  of  the  heat  transfer  for  the  region  near  the 
throat  which  is  clearly  incorrect. 

It  is  therefore  recommended  that  the  simple  Reynolds 
analogy  be  used  if  only  a rough  estimate  of  the  heat  transfer 
is  desired.  These  results  should  be  adequate  as  long  the 
wall  temperature  and  pressure  streamwise  gradients  are  not 
too  strong.  If  a more  refined  computation  is  necessary,  it 
is  felt  that  the  results  of  the  correction  procedure  (invoked 
by  HTCORR  = .TRUE.)  should  be  used  with  caution.  Also,  it  is 
helpful  to  remember  that  since  the  properties  of  the  velocity 
boundary  layer  depend  only  on  the  average  aspects  of  the 
thermal  layer,  the  computations  for  the  velocity  quantities 
may  be  accurate  even  though  those  for  the  heat  transfer  are 
not . 


E . Boundary  Layer  with  Transition 

A transitional  boundary  layer  case  consisting  of  flow 
over  a NACA  0012  airfoil  was  computed  as  shown  in  Fig.  21. 

This  example  was  taken  from  the  report  by  McNally  1 1 6 ] and 
was  used  in  explaining  the  input  procedure.  Section  III-D. 

The  data  and  predictions  for  the  displacement  and  momentum 
thicknesses  agree  very  well  except  for  just  downstream  of  the 
transition  point  where  the  computations  for  the  displacement 
thickness  show  a sudden  drop.  As  discussed  previously,  this 
is  a consequence  of  the  pointwise  nature  with  which  the  tran- 
sition process  is  assumed  to  occur  in  the  calculations.  Across 
the  transition  location,  H is  taken  as  continuous,  and  Z is 
chosen  to  make  the  momentum  thickness,  62,  continuous  (as  can 
be  verified  in  Fig.  21).  Other  boundary  layer  parameters  such 
as  displacement  thickness  and  skin  friction,  which  are  derived 
from  these  two  dependent  variables,  are  therefore  discontinu- 
ous at  transition.  The  comparison  of  the  data  and  computa- 
tions in  Fig.  21  for  the  displacement  thickness  suggests  that 
perhaps  a "hand-smoothing"  operation  is  warranted  to  remove 
the  discontinuity. 

Note  also  that  the  predicted  transition  location  agrees 
well  with  the  measured  transition  range. 


F . Miscellaneous  Examples 

The  final  three  examples  illustrate  the  way  in  which  pro- 
gram COMPBL  can  be  used  in  conjunction  with  inviscid  solutions 
to  predict  boundary  layer  characteristics  in  design  situations 
The  first  two  cases  consider  the  flow  over  the  10°  conical  af- 
terbody shown  in  Fig.  22.  The  inviscid  wall  Mach  number  dis- 


29 


tribution  was  determined  from  the  method  of  characteristics 
using  program  TSABPP-2  developed  by  Addy  [31].  The  approach 
Mach  number  is  assumed  to  be  = 2 and  the  ambient  static 
pressure  and  temperature  are  taken  as  sea-level  values, 

Poo  = 14.7  psia  and  Too  = 520°R,  which  yield  a very  high  unit 
Reynolds  number,  Reco/L  - 1.4  x 107  ft“l.  For  the  first  case, 
Fig.  22,  the  turbulent  boundary  layer  is  assumed  to  start  at 
the  expansion  giving  the  integral  parameters  shown  in  the 
bottom  half  of  the  figure.  Figure  23  shows  the  more  realis- 
tic situation  for  which  the  turbulent  boundary  layer  is  assumed 
to  begin  growing  on  the  missile  body  and  therefore  is  of  finite 
thickness  just  upstream  of  the  expansion.  This  second  case 
demonstrates  the  capability  of  program  COMPBL  to  calculate 
through  a "discontinuity"  in  the  Mach  number  distribution.  The 
wall  Mach  number  was  assumed  to  be  constant  at  M = 2 up  to 
X = -0.001  ft.  and  then  for  0 £ X £ 4 ft.,  the  method  of  char- 
acteristics results  were  used.  Note  that  the  expansion  fan 
produces  rather  large  discontinuities  in  the  shape  factor, 

6^/62 » and  the  skin  friction  coefficient,  Cf. 

The  last  example.  Fig.  24,  is  of  the  flow  through  a plane 
two-dimensional  Mach  2 nozzle.  The  wall  Mach  numbers  were 
computed  using  the  method  of  characteristics  nozzle  design 
program  NOZCS  written  by  Addy  with  typical  laboratory  operat- 
ing conditions  being  chosen  as:  PQ  = 600  kPa , T0  = 300  K, 
throat  radius  = 1 cm.  The  turbulent  boundary  layer  is  as- 
sumed to  start  at  the  throat.  Figure  5 shows  the  corrected 
boundary  coordinate  results  from  file  TAPE 4 . At  the  nozzle 
exit  the  corrected  centerline  distance  is  1.714  cm  compared 
to  the  method  of  characteristics  inviscid  design  value  of 
1.686  cm. 
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V.  CONCLUSIONS 


A FORTRAN  computer  program  based  on  Computational  Method 
II  of  Walz  |1]  has  been  developed  to  calculate  the  integral 
parameters  for  the  two-dimensional  boundary  layer  flow  of 
ideal  gases.  The  flexibility  of  the  program  allows  plane  or 
axisymmetric  geometries  with  the  boundary  layer  either  laminar 
or  turbulent,  incompressible  or  compressible,  and  either  with 
or  without  heat  transfer.  In  addition,  the  location  of  tran- 
sition can  either  be  predicted  or  specified  as  an  input. 

In  order  to  test  the  accuracy  and  reliability  of  the 
method  and  program,  a number  of  boundary  layers  were  computed 
for  which  either  data  or  analytical  solutions  are  available 
for  comparison.  Five  incompressible  turbulent  and  five  com- 
pressible turbulent  boundary  layer  examples  were  analyzed, 
most  of  which  fall  in  the  "difficult"  category.  For  most  of 
these  cases  the  agreement  between  the  data  and  the  predictions 
is  good.  For  the  laminar  and  transitional  boundary  layers  cal- 
culated, the  computed  results  agree  with  the  data  or  analyti- 
cal solutions  extremely  well.  However,  the  heat  transfer 
predictions  have  met  with  less  success.  For  rough  estimates 
of  the  heat  transfer,  it  is  suggested  that  the  default  Reynolds 
analogy  be  used.  The  heat  transfer  correction  procedure  or 
"modified  Reynolds  analogy",  discussed  in  Section  II-E,  should 
be  used  cautiously  since  it  has  been  only  moderately  success- 
ful . 

A major  asset  of  the  program  is  its  computational  speed. 

A typical  boundary  layer  can  be  calculated  on  the  CDC  CYBER 
175  in  \ - 1 CPU  seconds,  making  parametric  design  studies 
practical . 

Based  on  these  characteristics,  it  is  felt  that  program 
COMPBL  is  a useful  and  accurate  tool  for  estimating  the  inte- 
gral boundary  layer  parameters  in  a wide  variety  of  flow 
situations . 
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TABLE  I.  B vs  H from  Hartree  Solutions 


■ 


V 

B or  8e 

H 

-0.1988 

1.515 

-0.19 

1.5200 

-0.18 

1.5250 

-0.16 

1.5333 

-0.14 

1.5400 

-0.10 

1.5517 

0 

1.5720 

0.1 

1 . 5857 

0.2 

1.5950 

0.3 

1.6018 

0.4 

1.6070 

0.5 

1.6113 

0.6 

1.6150 

0.8 

1.6207 

1.0 

1.6250 

1.2 

1.6290 

1.6 

1.6345 

2.0 

1.6380 
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TABLE  II.  Summary  of  Procedure  for  Determining  Starting 
Values  ZSTART,  HSTART,  BSTART 


Sharp  Leading  Edge  or  Stagnation  Point 


ZS  TART : 


HSTART : 


plane  two-dimensional 

Sti  0 
3ti  = 0 

ax  isymivietr  ic 


0 < M <1 

oo  


M > 1 


BSTART : 


Table  I 


Fig . 2 


Table  I with  3E=jq3- 
Fig . 2 

(total  included  angle) 


(a)  AX1SYMMETRIC  CASE 


(b)  PLANE  TWO-DIMENSIONAL  CASE 


Figure  1.  Location  of  coordinates  X,  Y,  R;  approach  flow 
variables  (MINF)  , Re^/L  (REINFL)  , T^,; 

and  angle  S (BSTART) . 


i 
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Figure  2.  H vs  Mo>  used  in  determining  HSTART 
value  for  compressible  flows 
(after  Walz  11],  Fig.  6.5). 
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Figure  4(a).  Cone 


SPI.ACEMENT  AND  MOMENTUM 
THICKNESSES,  6.  and  6,  (ft) 


5,:  o McNally 

1 DATA 

PRESENT 

COMPUTATIONS 

6„:  A McNALLY 
* DATA 

PRESENT 

COMPUTATIONS 

IL  / 

TRANSITION 

RANGE 

So 

COMPUTED 

TRANSITION 

c/' 

jy'.s 

.A' 

i 

J L 

Figure  4(b).  Comparison  of  computed  results  and  data 
for  case  given  in  Figure  3. 


CORRECTED  BOUNDARY  CO-ORDINATES 


MACH  2 PLANE  2-D  NOZZLE  WITH  HYPERBOLIC  ENTRANCE  SECTION  (L^CN) 


AXIAL  LOCATION 

[L] 

I 


AXISYBBETPIC  B ADI  US,  OR 
2-D  DISTANCE  PROfl  CENTERLINE 

[L] 


• 

1. 00000 

8 12620L-01 

1 . 00099 

178443 

1. 00239 

264765 

1. 00525 

397739 

1.03948 

5 16029 

1.01522 

.0  199  1 

1.05423 

. 37ob3 

1. 09524 

. 7 101  3 

1.14248 

. 03971 

1 . 19671 

. 27076 

1. 23868 

. 00444 

1.36607 

. 79696 

1.48  018 

. 290  'j  b 

1. 54013 

. 6^512 

1. 59286 

.30462 

1. 63210 

. 60483 

1 . 66447 

. 26597 

1.66770 

. 76603 

1.70398 

'.  12769 

1.71067 

. 47646 

1.71 394 

Figure  5.  Example  printout  from  file  TAPE4 . 
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SKIN  FRICTION  COEFFICIENT, 


0.002 


0.001 


ol L 

100 


Figure 


J 1 

200 

X (in.) 


J L 

300 


8.  Comparison  of  computed  results  and  data 
for  Clauser,  flow  number  1 (2200). 


SKIN  FRICTION  COEFFICIENT,  C.  MOMENTUM  THICKNESS,  6,  (in.)  SHAPE  FACTOR,  S./8 
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X V 


Figure  12.  Comparison  of  computed  results  and  data 
for  Winter,  Smith,  and  Rotta  subsonic  waisted 
body  flow,  = 0.6. 
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SKIN  FRICTION  COEFFICIENT,  Cf  MOMENTUM  THICKNESS,  103  SHAPE  FACTOR,  6j/5 


Figure  13.  Comparison  of  computed  results  and  data 
for  Winter,  Smith,  and  Rotta  supersonic  waisted 
body  flow,  M„  = 1.4. 
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Figure  14.  Comparison  of  computed  results  and  data 
for  Winter,  Smith,  and  Rotta  supersonic  waisted 
body  flow,  M*,  = 2.8. 
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DISPI.ACEMENT  THICKNESS,  6j  (In.) 


PRESENT  COMPUTATIONS 

O ALSTATT  DATA  - PITOT 

A ALSTATT  DATA  - 
LASER  VELOCIMETER 


I 

0.3  f- 


X (in.) 


Figure  15.  Comparison  of  computed  results  and  data 
for  Alstatt-AEDC  transonic  flow. 
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a 


Figure  16.  Comparison  of  present  computations  and  finite 
difference  calculations  of  Smith  and  Clutter 
for  Hiemenz  cylinder  flow. 
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Figure  18.  Comparison  of  present  computations  and 
Flugge-Lotz-Eichelbrenner  finite  difference 
calculations  for  compressible  laminar  flow. 
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Figure  19.  Comparison  of  present  computations  and 
analytical  solution  of  Chapman  and  Rubesin. 
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X (cm) 


Figure  20.  Comparison  of  present  computations  and 
data  of  Boldman,  Schmidt,  and  Ehlers. 
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Figure  21.  Comparison  of  present  computations  and  data 
taken  from  McNally  [16]  for  flow  over  a NACA  0012 
airfoil.  (This  is  a repeat  of  Figure  4(b)). 
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PROGRAM  COHPBL  (INPUT,  OUTPUT,  TAP  E4,  TAP  £5- IN  PUT,TAPEb=OUTPUT) 

WRITTEN  BY: 

J.C.  DUTTON 

DEPARTMENT  OF  MECHANICAL  AND  INDUSTRIAL  ENGINEERING 
UNIVERSITY  OF  ILLINOIS  AT  URBAN A- CHAMPAIGN 
URBANA,  ILLINOIS  61801 
OCTOBER,  1877 

PROGRAM  COMPBL  IS  A TWO-DIMENSIONAL  BOUNDARY  LAYER  ANALYSIS  PROGRAM 
FOR  IDEAL  GASES  BASED  ON  COMPUTATIONAL  METHOD  II  OF  WALE  ("BOUNDARY 
LAYERS  OF  FLOW  AND  TEMPERATURE”,  MIT  PRESS,  1969)  AND  THE  DISSIPA- 
TION LAWS  OF  FELSCli  AND  ROTTA-TRUCKENBRODT.  THE  METHOD  HAY  BE 
APPLIED  TO  INCOMPRESSIBLE  AND  COMPRESSIBLE,  LAMINAR  AND  TURBULENT, 
PLANE  2-D  AND  A XISYMMETRIC  BOUNDARY  LAYERS,  WITH  AND  WITHOUT  HEAT 
TRANSFER.  THE  COMPUTATIONS  MAY  BE  STARTED  AT  AN  ARBITRARY  STREAM- 
WISE  LOCATION  WITH  THE  BOUNDARY  LAYER  EITHER  LAMINAR  OR  TURBULENT, 
AND  TRANSITION  AND  SEPARATION  ARE  PREDICTED.  THE  INTEGRAL  MOMENTUM 
AND  MECHANICAL  ENERGY  EQUATIONS  ARE  SOLVED  ITERATIVELY  IN  A STEP-BY- 
STEP  FASHION  WITH  THE  FREESTRSAfl  VELOCITY  APPROXIMATED  WITH  A 
PIECEWISE  LINEAR  FUNCTION.  FOR  CASES  WITH  HEAT  TRANSFER  AN  OPTION 
IS  AVAILABLE  (UTCORR=. T RU E. ) WHEREBY  THE  CALCULATION  OF  THE  WALL 
HEAT  FLUX  IS  IMPROVED  BY  SOLVING  THE  THERMAL  ENEBGY  INTEGRAL 
EQUATION  FOR  A CORRECTION  PARAMETER,  CUT,  TO  ACCOUNT  FOR  THE  EFFECTS 
OF  STREAMWiSE  WALL  TEMPERATURE  ANE  PRESSURE  GRADIENTS. 

THE  DIMENSIONAL  VARIABLES  IN  THE  ANALYSIS  ABE  Z START, REI NFL , XTR AN , 
XI,YI,FI,ZI,D1,D2,B999,  AND  RCORR  WHICH  ALL  HAVE  DIMENSIONS  OF 
LENGTH,  [L]OR  (L]**(-1).  iS TA RT  , REI NFL, XTR AN, XI , YI , AND  RI  MUST 
BE  INPUT  WITH  CONSISTENT  UNITS,  AND  THEN  ON  OUTPUT,  ZI , D 1 , D2, D9 99, 
AND  RCORR  WILL  HAVL  THE  SAME  UNITS. 

THE  INPUT  VARIABLES  FROM  FILE  INPUT  ARE: 

FLOW LITERAL  VARIABLE  DESCRIBING  FLOW  REGIME — EvUAL  TO 

"LAMINAR”  Oh  "TURBULENT" 

GEOM LITERAL  VARIABLE  DESCRIBING  GEOMETRY — EQUAL  TO 

"PLANE  2-D"  OR  "AXISYM" 

MTYPE LITERAL  VARIABLE  DESCRIBING  VELOCITY  INPUT  DATA  TYPE 

— EQUAL  TO  "MACH"  OR  "MSTAR" 

TTYPE LITERAL  VARIABLE  DESCRIBING  WALL  TEMPERATURE  INPUT  DATA 

TYPE — EQUAL  TO  "THETA"  CR  "TWTI NF" 

ZSTART-- STAR TING  VALUE  OF  Z,  JL]  (DEFAULTED. D) 

HSTAHT — STARTING  VALUE  OF  H;  SEF  WALZ  PP  2b7-268  AND  FIGS  1.1 
AND  6.5  AND  ACCOMPANYING  REPORT,  TABLES  I AND  II 
AND  FIG.  2 (DEFAULT1  1.  57  2) 

BSTART — BSTAhT*PI  IS  THE  INCLUDED  ANGLE  AT  THE  LEADING  EDGE 
FOR  BOTH  PLANE  2-D  AND  AXISYMMETR1C  EXTERNAL  FLOWS 
(DEFAULTED. 0) 

MINP MACH  NUMBER  (OR  MSTAR)  OF  APPROACH  FLOW 
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10 

COM 

20 

COM 

30 
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40 
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53 
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60 
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70 
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80 
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90 
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COM 

110 
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120 

COM 
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140 
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160 
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190 

COH 

200 
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210 
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230 
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24  3 

COM 
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COH 
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COM 

280 

COM 
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COM 

300 

COM 

310 

COH 

320 

COM 

330 

COM 

340 

COH 

353 

COM 

36  3 

COM 

370 

COH 

38  3 

COH 

390 

COM 

430 

COM 

410 

COM 

423 

COM 
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COM 

443 

COH 

450 
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COM 
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COH 
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PBOGBAH  COHPBL  . 


C...  BEI NPL — 8 EYNOLDS  NUMBER  DIVIDED  BY  C H ABACTEBI STIC  LENGTB  OP 
C...  APPBOACH  FLOW,  I.  8.  BH01NP*UI NP/HUINF, 

C...  XTBAN X LOCATION  POR  SPECIFIED  TBANSITION,  £L] 

C...  (DEFAULT-0.0) 

C...  PSTINT — FREESTREAM  TURBULENCE  INTENSITY — USED  IN  TBANSITION 
C...  SUBROUTINE,  IN  PERCENT  (DEPAULT-O.  3) 

C...  G RATIO  OF  SPECIFIC  HEATS  (DEFAULT =1.405) 

C...  B EXPONENT  ON  VISCOSITY  POBEB  LAB  (DEPAULT-O.  7) 

C...  RL LANINAB  RECOVERY  FACTOR  (DEF AULT-0. 85) 

C...  RT TURBULENT  RECOVERY  FACTOR  (DEF AULT-0. d8) 

C...  SL LANINAB  REYNOLDS  ANALOGY  FACTOR  { DEFA ULT-0. 80) 

C...  ST TURBULENT  REYNOLDS  ANALOGY  FACTOR  (DEFAULT-0.82) 

C...  EPS CONVERGENCE  CRITERION  VARIABLE  (DEFAULT- 1. OE-4) 

C...  NOTRAN — LOGICAL  VARIABLE  WHICH  IF  .TRUE.  SUPPRESSES 
C...  CALLING  OF  THE  TBANSITION  SUBROUTINE  FOR  LANINAB 

C...  BOUNDARY  LAYERS  (DEF A ULT- . FA LSE. ) 

C...  HTCORR — LOGICAL  VARIABLE  WHICH  IF  .TRUE.  CAUSES  THE  THERNAL 
C...  ENERGY  INTEGRAL  EQUATION  TO  BE  SOLVED  FOR  THE  HEAT 

C...  PLUX  CORRECTION  PARANETEB,  CUT  (D EFAULT-.F ALS E.) 

C...  EBSCR LOGICAL  VARIABLE  WHICH  IP  .TRUE.  CAUSES  INTERN  EDI ATE 

C...  H VALUES  AND  VARIABLES  ASSOCIATED  WITH  THE  TURBULENT 

C...  DISSIPATION  INTEGRAL  AND  HEAT  TRANSFER  CORRECTION 

C...  PARAMETER  TO  BE  PRINTED  FOB  DEBUGGING  PURPOSES 

C...  (DEFAULT-.  FALSE.) 

C...  XI AXIAL  LOCATION,  [L] 

C...  Y I NORMAL  LOCATION,  [L] 

C...  NDI FREtSTREAN  MACH  NUMBER  (OR  NSTAR) 

C...  RI CROSS-SECTIONAL  BADIUS  OF  AXISYNNETBIC  BODIES,  OR 

C...  LOCATION  NORMAL  TO  CENTERLINE  FOR  2- D BODIES,  [L] 

C...  1HETI WALL  TEMPERATURE  RATIO,  (TADW ALL-T WALL) / (TA DBA LL-TSTBEAN) 

C...  (OR  TWALL/TINF) 

C 


C...THFSE  VARIABLES  SdOULD  BE  INPUT  IN  THE  FOLLOWING  WAY.  THE  FIRST 
C. . . CARD  (RECORD)  IS  FOB  A TITLE  TO  HELP  IDENTIFY  THE  OUTPUT.  ANY 
C... MESSAGE  UP  TO  80  COLUMNS  CAN  BE  USED,  BUT  FOR  AESTHETIC  REASONS 
C ...  THE  MESSAGE  SHOULD  BE  CENTERED  IN  THE  80  COLUMNS.  IT  IS  ALSO 
C. . . SUGGESTED  THAT  THE  UNITS  FOR  DIMENSION  [L]  BE  ENTERED  AS  PART  OF 
C ...  THE  TITLE.  ON  THE  NEXT  CARD  THE  FOUR  LITERAL  VARIABLES  FLOB.GEON, 
C...MTYPE.TTYPE  ARE  ENTERED  IN  «A 1 0 FORMAT.  THE  VALUE  OF  EACH  OF 
C... THESE  LITERALS  SHOULD  BE  LEFT  JUSTIFIED  IN  EACH  10  COLUMN  BLOCK. 

C. . . NAMELI ST  "DL"  WHICH  ENCOMPASSES  VARIABLES  2START-ERBOH  IS  ENTERED 
C...ON  THE  NEXT  CAHD(S).  THE  FIRST  COLUMN  ON  THE  NAMELIST  CARD  (S) 

C... SHOULD  BE  BLANK.  ALL  VARIABLES  IN  NAMELIST  "BL"  EXCEPT  HINF  AND 
C... RE  IN  F L AHE  DEFAULTED,  SO  ONLY  THESE  TWO  AND  THOSE  WHICH  ARE  TO  BE 
C... OVERRIDDEN  NEED  BE  ENTERED.  THE  REMAINING  CARDS  CONTAIN  THE  LOCAL 
0 . . . DATA : X.Y.M.R, THETA  FOB  EACH  POINT  AT  WHICH  THE  BOUNDARY  LAYER 
C... PARAMETERS  ARE  TO  BE  FOUND.  THESE  VARIABLES  ARE  ENTERED  IN  5F10 
C... FORMAT,  AND  AS  MANY  POINTS  AS  DESIRED  CAN  BE  USED. 

C 
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PBOGBAN  CONPBL  . 


C...THE  OUTPUT  VARIABLES  TO  PILB  OUTPUT  ABE: 


C 

C...  XI AXIAL  LOCATION,  [ L] 

C...  ZI ZI»D2« <RD2**N)  (N*1.0  FOB  LANINAS  BOUNDARY  LAYERS; 

C...  0.26b  FOB  TUBBULENT  BOUNDABT  LAYEBS),  [L] 

C...  HI HI=  (D3/D2)  U (D  3s  BN  ERG  Y LOSS  THICKNESS) 

C...  D1D2 SHAPE  FACTOR,  D1/D2 

C...  D 1 DISPLACEMENT  THICKNESS,  [L] 

C...  D2 HONENTUN  THICKNESS,  [L] 

C...  D999 99.9*  BOUNDARY  LAYER  THICKNESS, [ L ] 

C...  RD2 SOHENTDH  THICKNESS  REYNOLDS  NUMBER  (BHOD*UD*D2/NUHALL) 

C...  CF LOCAL  SKIN  FRICTION  COEFFICIENT  (2*TA  UH  ALL/ (RHOD*  UD**2)  ) 

C...  QDIH DIMENSIONLESS  LOCAL  HALL  HEAT  FLUX  (QH ALL/ (RHOD*U D** 3 ) * 

C...  -R*STANTON*THETI/2) 

C 

C...THE  OUTPUT  VARIABLES  TO  FILE  TAPE4  ABE: 

C 

C...  XI AXIAL  LOCATION,  £ L ] 

C...  RCORR CORRECTED  RADIUS  OR  NORMAL  LOCATION — EQUAL  TO  RI*D1,  [LJ 

C 


IMPLICIT  REAL  (J,R,N) 

LOGICAL  ERROR, NOTRAN, HTCOBR 

COMMCN/CTRANS/RD2,  N,  H , RT,  S,  ST,  FLOW.G,  HTYPE,  TTY  PE,  MI  NF,  KTRANS, 
$KINST,RD2U,PPARI,PPABIM1,DE LX, F STINT 
01  MENSICN  HMAIN  (21)  , HINNEB  (21)  , QHKITE  (8) 

...DECLARE  INPUT  NAMELIST  AND  SET  DEFAULT  VALUES 

NA  ME  LI ST/BL/  ZSTART. HSTART, BSTA RI , MIN F , REINFL, ITR AN, FSTINT, G, R , 
SPL ,RT, SL, ST, EPS, ERROR, NOTHIN, HTCOBR 
DATA  ZSTART, HSTART, BSTART.CSTART, XTRAN , FSTINT , G, H, R L, RT.SL,ST,EPS 
i ERROR, NOTRAN, HTCORR/0. 0, 1. 572, 0.0, 0.0 ,0.0, 0.0, 1.405,0. 7,0.85, 

$0. 88,0.80,0.82, 1.0E-4, .FALSE. , . FA LSE. ,. FALSE./ 


...GENERAL  STATEMENT  FUNCTIONS 

FMS  (G,M)  =SQRT  l ( G ♦ 1 - ) /2 . *M**  2/  (1  .♦  (G- 1 . ) /2.  * M**  2)  ) 

FTNT1NF  (R,G,  MD,  MI  NF.  THETA)  = ( 1 . ♦ R*  (G-  1 . ) /2.  *NU»»  2*  ( 1 .- THET  A)  ) * 
i (1 . ♦ (G-  1.)/2.*MINF**2)/(I.*  (G-  1 . ) /2.  * MD**2) 

FB  (THETA,  R,G,  HD)  =THETA*R*  (G-  1.)  /2.  *MD**2 

FC1  (REINFL,  MDS,  MI  NFS,  MD,  MIN P,R,  G,  TiitT  A,  H)  -REINFL*  MDS/MINFS»  ( ( 1 . ♦ 
i (G-1.J/2.  *NINF**2)/(1.+  (G-1  .)  /2.*MD**2)  ) «*  (1./(G-1.)  ) * ( (!.♦  (G-1.  ) 
$/2.  *MD»*2)/(  ( 1 . *R*  (G-  1.  )/2.  *M  D*  *2«  ( 1 . -TH  ET  A)  ) * ( 1 . ♦ (G-  1 . ) /2  . *MINF 
$»*  2)  ) ) «♦« 

FD  2 (N,Z,C1)  = (Z/C1**N)  **  ( 1 ./  ( 1 . ♦ N)  ) 

FR  02  (N,Z,C1)  = (Z*Cl)»»  (1./(1.*N)  ) 

FRD2U  (RD2,D2D2U,R,  G,MD,  THETA)  =RD2/  (D2  D2 U*  ( 1 . *R»  (G- 1 . ) /2 . *MD**2* 

$ (1  .-THETA)  ) ) 

FCF  (A,RD2,  N,  D2D2U)  = 2. * A/RD2**  N*  D2  D2U 
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COM  993 
COM 1000 
COH 1 010 
COH1020 
COH 1030 
COH  104  0 
COH 1050 
COM1060 
COM1 070 
COM  108  0 
COH 1090 
COH1 100 
COH1 110 
COH1 120 

coni  no 

C0M1 140 
com  150 
C0M1 160 
COM  1 170 
COM  1 18  0 
COH1 190 
COH  1200 
COM  1210 
COM1 220 
COH1230 
COM  1240 
COH 1250 
COM1 26  0 
COM1270 
COM  1 28  0 
, COM  129  3 
COM  1300 
COH 1310 
COH1 320 
COM  1 33  0 
COH1340 
COM  1350 
COM1 360 
COM  1370 
COM1 380 
COH1 390 
COM  143  0 
COH 1410 
COM  1420 
COH 1430 
COH 1440 
COM  145  0 
COH 1460 
COH 1 47  0 
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FQDIN  (CP,S,B,CHT,G,MD)=-CF/(2.*S)*  (B*CHT)  / (HD**2*  (G-1.)  ) 

PDELX  (XI,  YI.XIMI.YIMI)  = SQRT  ( ( XI-X  IM  1)  **2*  (II- XIM1  ) • *2) 

FPPAH {C1,DNSDXHS,D20,R,G,HD,THETA)=C1*D8SDXHS*D20**2/(1.*a*{G-1.)/ 
*2.  *HC**2*  (1. -THETA) ) 

...UNIVERSAL  FUNCTIONS  FOR  THE  METHOD 

PAZ  (UIM1UI.F1B,  NUMRN)  = UIN1UI**F1B*MUtiRN 

FBZ  (UIM1UI,P1B,MUNBRN)= (1.-UIH1UI**  ( 1 . *F  1 B)  ) / ( ( 1.  ♦ F IB)  * (1  .- UIM 1UI 
S)  ) *NUHBRN 

FAH  (UIM  1UI  ,F3B)  =UIN1UI**FJB 

PBli  (UIM  1UI,  P3B)  * ( 1.  — UIM  1UI*  * ( 1.  *F  3B)  ) / ( ( 1.  *F3B)  * (1.  - UIM  101)  ) 

FF1 (N,D1D2,MD)=2.*N*(1.*N) *D1D2-MD**2 

FF2  (N,  D2D2U  r A)  = ( 1 . ♦ N)  * D 2D2U*A 

FP3(D1D2,  D4D3)  = 1 . -D  1D2*2. *D4D3 

PP4 (RD2,N,CD,A,D2D2U#HS) =BD2**N*CD- A* D2 D2U*HS 

FF5  (DNSDXNS,hS,D1D2,G,MD,HSP,  DD  2DXD  2)  = (DMSDXH  S*  (1.*D1D2* 

S (G-J.)  *MD**2*HS*(1.-G*HD**2))  *HSP*HS»  DD2DXD2)  / (jiS-  1 . ) 

FF6 (HS,BP,B,DMSDXMS,D1D2,G,ND)= (BP-B*DM3DXMS* (1.*D1D2* (G-1.) * 
$ND**2))/ (HS-1.) 

FD1D2  IH  12,  D2D2U,R,G,ND,  HS, THETA  ,CHT)  = H 1 2/D2D2D*R*  (G-1  . ) /2.  *HD**2  * 

$ (HS-THETA)  *CHT*  (HS-1.) 

FD4D3 (R,G,  MD  , HS,  TH  ETA,  CHT)  = H*  (G-  1 . ) /2.  » MD** 2*  (HS-TH  BT  A) /HS 
f*CHT*(HS-1.)/HS 

FD2D2U  (R,G,ND,HS,TilSTA)  = 1./  ( 1 . ♦ B*  (G- 1 . ) /2.  * HD*  • 2*  (HS-THETA)* 

$ (2.-HS)  ) 

FHS(H,  PSI)  =H*PSI 
Fll  (H  S,  PSI)  =HS/PS1 

FPSI  (PSI12,MD,FS1JP)  * 1.  ♦ (PSI1  2—  1.  ) *MD/  (MD*  (PS1 1 2- 1.  ) /PSI  JP) 

FPSI  12  (D10D,H,THETA,G1)  = (2.-D1UD)  *THETA/H*  (1.-D1UD)  * (1.- 
STHETA)/  (11*G1) 

FPSIOP  (H,  THETA)  =0. 0144*  (2.-H)  * (2. -THETA)  **0.8 

...STATEMENT  FUNCTIONS  FOR  LAMINAR  BOUNDARY  LAYERS 
FAL(H)  = 1.72b1*  (H- 1 . SI  5)  *•  3.  71  58 

FCDL  (H,R,G,HD,  THETA,  H,  RD2,  D2D2U)  = 2.  *D  2D2U/BD2*  (0.  1564*2.  1921* 

$ (H- 1.5  15)  **1.  73)  * ( (1.  *R*(G-  1.)/2.  *MD**2*  ( ( 1 . 1 6 3*H- 1 . 372) -THET  A* 

* (2.*d-2.581) ) )/(1.*R*(G-1.)/2.*HD**2* (1. -THETA) )) **W 
?H  121(H)  =4.0306-4.284  5*  (H- 1 . S 15)  **  3.  388b 
FD1UDL  (II)  = 0.  4 20-  (H-l.  515)  **  (0.  424*H) 

FG1L  (H)  =0.324*0.336*  (H-l.  51  5)  **0.555 

...STATEMENT  FUNCTIONS  FOR  TURBULENT  HOUNDARY  LAYERS 
FAT  (H)  =0.0389  4*  (H- 1.  5 1 5)  *«0 . 7 

FCDTF  (CP,RD2,N,JP,J,JSP,HS,D1D2,PI)  =CP/2.*  ( (J  P-N/2.  *J)  / (R  D2** 

J (1/2. ) *JSP)  * ( 1.*(D1D2*1.)  *PI/D1D 2)  *HS*  (1.  ♦ (D1D2-1.)  /D1  D2*PI) ) 
FCDTRT  ( J2D2U , RD 2)  =0.0112*  D2D2  U*RD2**  (-0.  168) 
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COM  1 490 
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COM  1510 
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CON1540 
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APPENDIX  A.  COMPRESSIBLE  BOUNDARY  LAYER  PBOG BAN. • . COHPBL 

PROGRAM  CONPBL  . 


FJ  (D1D2,CF)=  (D1D2-1.)/(D1D2*SQHT  (CP/2.)  ) 

FA  1 (D1D2|  = 0.  029  *(0.93 -1.95*  ALOG  1 0 ( D 1 D2)  ) **1.  705 

FJSP  (A  1,D1D2,HS)  = (D1D2-  (D 1D2-  1.  ) ♦ ( 1.- 0.  0209366/A  1*  (0.  93-1.95* 

$ ALOG  10  (D1D2) ) **0.7 05)  ) / (SQRT  ( A 1 ) *D1D2**2)  * (0. 5521 6-0. 3875*HS* 

*0.  04855*SQBT  (0.  775*HS-  1.  10667)  ) / ( (US-  1.  43  1)  **2*SQRT  (0.  775*HS- 
$1.  10667)) 

FPI  (D1D2,D2,CF,  UHSDXflS)  =- 2.  *D  1D2*D2*DHSDXHS/CP 
F I1 1 2 T { 8 ) - 1 . ♦ 1 • 48*  (2.-H)  *104.*  (2.-H)  **6.7 
FD1UDT  (U)  = (2.-11)/  (2.  *(H-1.)  ) 

FG  IT  (8)  =0.  30  6*  (U- 1.5)  -0.885*  (H-  1.  5)  **  1.  53 

C 

C...READ,  CHECK,  AND  WhITE  INPUT  VARIABLES  AND  SET  INITIAL 
C.. .PARANETEES 

C 

IP=KTRANS=KINST=LCOUNT=LIN=0 

RRA  T 1PN=  1.  0 S DEL  X= DELI IN  1=  0.  0 

READ  (5,900)  Q WRITE  , FLOW  ,G  EOH,  NT  YPE,  TT  YPE 

RE  AD  (5 , BL) 

K=  0 $ WRITE  (6, 901)  $ WRIT£(4,901) 

IF  (FLOW  . EQ.  "LAMINAR"  .OR.  FLOW  .EQ.  "TURBULENT")  GO  TO  10 
K=  K* 1 i W RITE (6,902) 

10  IF (G EON  .EQ.  "PLANE  2-D"  .08.  GEOM  . EQ.  "AXISYS")  GO  TO  20 
K=  K*  1 $ W RITE  (6,903) 

20  IF  (MTYPE  .EQ.  "MACH"  .OR.  (1TY  PE  . EQ.  "MSTAR")  GO  TO  30 
K=  K* 1 $ W RITE (6,904) 

30  IF (TTY  PE  .EQ.  "THETA"  .OR.  TTYPE  .EQ.  "TWTINP")  GO  TO  40 
K=  K*  1 S WRITE (6 , 905) 

40  I P (USTART  .GT.  1.515  .AND.  HSTART  . LE.  2.0)  GO  TO  50 
K=  K ♦ 1 $ W RITE (6,906) 

50  IF  (BSTART  .GE.  -0.  1988  .AND.  B SI  ART  . LE.  2.0)  GO  TO  60 
K=K+1  $ WRITE  (6,907) 

60  IFIZSTA8T  .GE.  3.0  .AND.  MINE  . GT.  0.0  .AND.  REINFL  .GT.  0.0 

S.AND.  FSTINT  .GE.  0.0  .AND.  W . GE.  0.0  .AND.  RL  -GT.  0.0  .AND.  R 
$.a T.  0.)  .AND.  SL  .GT.  3.0  .AND.  ST  .GT.  0.0  .AND.  EPS  .GT.  0.0) 
$GO  TC  70 

K=  K»  1 $ W RI TL  (6 , 9 38) 

70  IF (G  .GE.  1.0  .AND.  G .LE.  1.67)  GO  TO  80 
K=KH  S W RI  TE  (6,909) 

80  IF  (.  NOT.  (MTYPE  .EQ.  ".1STAR"  .AKC.  NINP  .GT.  SQRT  ( (G  + 1.  ) /(G-  1.  ) ) ) ) 
$GO  TC  90 

K=K*1  i WRITE  (6,9 10) 

90  IF  (K  .EQ.  0)  GO  TO  100 

WRITE  (6 ,91  1)  i GO  TO  600 
130  CALL  MSOKT  (MINF,  MINFS) 

WRITE (6, 912)  Q W RI T E, P LO W, GEOM, B T YP E, TT Y PE, ZST ART, HSTART, BSTART, 
$NINF, REINFL, XT  RAN, FSTINT, G , W, RL,KT, SL.ST, EPS, NOTR AN ,HTCOR R, ERROR 
WRITE  (4,913)  Q WRITE 

IF  (FLOW  .EQ.  "TURBULENT"  .AND.  ZSTART  . EQ.  0.0)  KTR ANS= 1 
IF (FLOW  .EQ.  "TURBULENT"  .AND.  KT  HANS  .NE.  1)  GO  TO  110 
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APPFNDIX  A • COBPB  ESSIBLE  BOU  N D AB  i LAYER  PROG  RAB...COBPBL 

PROGRAM  COB  PB  L . 


N=  1 . 0 * a=RL  $ S=SL 

FLOM="LABINAB"  i GO  TO  120 

119  N=  J. 268  $ B = RI  S S»  ST 

...STARTING  SEQUENCE  FOR  BOUNDARY  LAYER  CALCULATIONS 

120  WRITE  (fa. 914) 

WRITE(4,915) 

IF (ZSTART  .NE.  0.0)  GO  TO  150 

IF  (G  ECB  . Eg.  "AXISYB")  BSTA RT=8STAR T/ ( 3. - BSTABT) 

NST A RT=  BSTART/  (2.-BSTART) 

READ (5,916)  XI, YI.BDi, Hl.THETI 
CALL  J1SCRT  (BDI.BDIS) 

RCOHR=RI  $ IP=IP+1 

IP  (SCI  .EQ.  0.0)  GO  TO  130 

D1=D2=D999=RC2= 3. 3 J DI D2=CF=QD1B=" " 

WRITE  (6,  917)  XI , Z ST  ART,  HSTA  RT  , D 1D2.D1  , D2 , D999 . R D2 ,CF, QDIS 
GO  TC  140 

130  RD  2=  0.  0 $ D1=D2=D999=D1 D2=CF=qDIM=" " 

WRITE (6, 918) XI.ZSTART.HSTABT, 01 02, D1, D2 , 0999, R D2 , CF , QOI B 
140  WRITE  (4,919)  XI.RCORR 

XI M 1=  XI  $ YIM1=YI  i B DI B 1 = SDI 

BDIB1S=ilDIS  i PPARI  B 1 = 0 . 0 

150  FEAD (5,9  16)  X 1 , YI , B DI , R I , TH  ETI 
CALL  BSORT  (BDI.BDIS) 

CALL  TSORT  (THETI.TWTII.BI.MDI)  i KPOI NT= 1 

IF  (Z  START  .NE.  3.3)  GO  TO  163  $ D ELX=  FDELX  ( XI  , Y I,  XI B1  , YIB  1 ) 

DSSDXBS*  (BDIS-MDIM  IS)  / (DELX*FBS  (G  , (BDI*  BDIB 1) /2.  ) ) 

CALL  1ERROR  <D  ELX  , SD1 , BDIS  ,G  , RI , GLOB  , T WT1I)  , RETURNS  (60  3) 

160  ilI  = H START  I ZI=ZSTART  $ CUT=CSTART 

IF  (FLOW  . LQ.  "TURBULENT")  GO  TO  17) 

A=  FA  L (ill)  t H12=FH12L(HI) 

D1UD=FD1UDL(liI)  i G1=FG1L(HI)  $ GO  TO  1 BO 
170  A=FAT(rii)  i li  1 2=FU  1 2T  (H I) 

DI  UD  = FD  1 UDT  ( HI)  $ G1=PG1T(H1) 

180  IF (ABS(THETI)  .LE.  2.0)  GO  TO  190 
PS1= 1.3  i GO  TO  230 
190  PSI0F=FPS10P  (Hl.THETI) 

PS  I 1 2=  FP  S1 1 2 (DI  UL, Hl.THETI,  G1) 

PS  1=  FP  Si  (PSi  12.BDI,  PSIOP) 

200  HIS=FHS  (HI, PSI) 

D2D2U  = FD2D2U  (R.G.BDI,  HIS.TH  ETI) 

DI  D2  = FD  1D2  (H  1 2,  D2D2U  , K , G,  MDI,  HIS,  THET1 , CHT) 

IF  (ZSTART  .EQ.  0.0)  ZI=  PF  2 ( N,  D2D2U,  A)  *DEL  X/  ( 1 . +BSTART*  FF  1 (N  , D 1 D2 
$BDI)  ) 

IF (Z  START  .EQ.  3.3  . AND.  GEON  . EQ.  "AXISYB")  ZI  = ZI/3. 

C 1 = FC  1 (REINFL.BDIS,  BINPS,  BDI,  BINE,  R.G.THLTI,  W) 

GO  TC  560 
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...READ  AMD  CHECK  INPUT  DATA  POB  MEET  POINT 

210  READ  (5,9  16)  XI,  Yl,  MDI , RI , TH  ETI 
IF  (EOF  (5))  600,  220 
220  CALL  MSOHT  (MDI,  MDIS) 

CALL  TSORT  (THETI,  THTII , Bl  , MDI) 

DELX=PDELX  (X1,YI,XIM1,YIM1) 

CALL  I ERROR (DELX, MDI, MD IS ,G , RI , GEOM, TMTI1)  .RETURNS  (600) 
KPCINT=KPOINT*1  S KTRIP=0  $ GO  TO  250 

...LINEAR  INTERPOLATION  FOR  SUBDIVIDED  POiNTS 

230  Di  V=  FLOAT  (LIM) 

DX?  = (XI-XIM1J/DIV  J DY1=  (YI-YIM1J/DIV 

DM  1=  (MDI-HDIM1)/CIV 

DT1=  (THETI-THETIM  1) /DIV  $ UR  1 = (HI-RIMl)  /DI  V 

260  XI=XIM1*DX1  f YI=YIM1*DY1  $ MDI=MDIM 1 *DM 1 

THETI^THETI  M 1 + DT 1 $ HI=RIN1+DR1  $ MDIS=FMS  (G,  MDI) 

BI  = FB  (THETI, R,G, MDI)  f TUT  1 1=  FTHTINF  (H  , G , MDI , MI  NF , THETI) 

IF  (FLOW  . EQ.  FLOWIM1  .OH.  ICOUNT  .hi).  1)  GO  TO  260 

...ALLOW  FOR  INTERVAL  SUBDIVISION  IF  NECESSARY 

250  ICOUNT=LCOUNT=LIM=0 
260  ICCUNT=ICOUNT+ J 

UIM1U1=MDIM1S/MDIS 

IF  (ICOUNT  .El).  1)  LIM  1 = IFIX  (AES  (1./UIN1UI-  1.)  /0.  0 15)  *1 

...CALCULATE  SOME  CUANTITIES  WHICH  ARE  CONSTANT  THROUGHOUT  THE 
...ITERATIONS 

DEIX-FDELX  (XI,YI,XIH1,YIH1)  t MDB=  ( MDI  *MDIN  1)  /2  . 

MDSB  = FMS  (G.MDB)  i TH ETB  = ( TH  El  i*THETI M 1 ) /2. 

rWTIB=FTWTiNF(R,G,MDB,MINF,THETE)  $ BB=FB(TH£TB,R,G,KDB) 
DM  SDXMS=  (MDIS- MDI  MIS)  / (DELX*MDSB) 

BP=  (3I-BIM1)  /DELX 

MU WR N-  (TWTIIM  1/TKTII)  **  (W*N) 

M U W B R N = (TWTIB/TWTi  I ) **  (W*M) 

IF  (G  EOM  -El).  "AXISYM")  RR  AT  1 PN=  ( RI/Iil  Ml ) * * ( 1.  ♦ N) 

C1B=FC1  (REINFL,  MDSB,  MI  NFS  , M DB , M IN  P,  R,  G,  THETB,  W) 

C 1 = FC1  (REIN PL, MDIS, MI  NFS, MDI, MI NF , R , G , THETI , W ) 

...MAIN  ITERATION  LOOT  FOR  BOUNDARY  LAYER  CALCULATIONS 

1=  KCDS  = KTRIP=0  i DilDX  = DHDXl  H 1 $ CHT=CHTIM1 

IF  (D  ELXI M 1 .EQ.  0.0)  GO  TO  270 
IF  (ABS  (DELX/DELXIM1)  .GT.  2.0)  CUDX=0.0 
270  HI0LD=HMA1N (1)=HIM 1 ♦DHDX*DELX 
IF  (ERROR)  WRITE  (6 , 920)  HMAIH(I) 
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PROG  BAfl  CONPBL  . 


280  1=1*1  $ HB  = (U10LD*HIN1)  /2  • $ CHTB=  (CHT*CHTIN1) /2. 

rE(EBBOB)  BBJTE  (6,921)  88 
IP  (HB  .LE.  1.515)  GO  TO  510 

IP  (HB  .LE.  2.0)  GO  TO  290  $ W BITE  (6, 922)  XI  $ GO  TO  600 

290  IF  (FLOW  .EQ.  "TURBULENT")  GO  TO  300 

...EVALUATE  AVERAGE  QUANTITIES  FOB  LANINAS  PLOWS 

AB  = P AL  ( H B)  $ II12B=PH12L  (HB) 

D1  UDB=FB1UDL  (HB)  $ G 1B=  PG  1 L (HB)  $ GO  TO  310 

...EVALUATE  AVERAGE  QUANTITIES  POB  TURBULENT  PLOWS 

300  !_  = F AT  (HB)  S H 1 2B=  PH  1 2T  (H  B) 

D1UDB=FD1UDT  (HB)  $ G 1B=  PG  IT  (H  B) 

...EVALUATE  UNIVERSAL  FUNCTIONS  IN  OBDEB  TO  OBTAIN  A NEW 
...VALUE  FOR  Z AT  STATION  I 

313  IP  (AES(THSTB)  . LE.  2.3)  GO  TO  320 
PSIB= 1. 0 $ GO  TO  330 

323  PSI3PB=FPSI0P(HB, THETB) 

PS1 1 2B=  FPS1 1 a (D1UDB,HB,THETB,G1B) 

P3IB=PPSI (PSI12B,NDB,PSI3PB) 

330  HSE=FHS (HB.PSIB) 

D2D2UB=FD2D2U ( B, G , NDB, HSB .THETB) 

04  C3B=FU4D3  (B,«,HDB,  HSB  , THETB  ,CHTB) 

DlD2fi=f D1D2(H12B, D2D20 B , H , G , NDB , HSB , THETB, CHTB) 

F1B=FF1  (N,D1D2B,NDB)  t F2B=  FF2  (N  , D2D2UB,  AB) 

P3B=FF3 (D1D2B.D4D3B) 

BZ=HUWBRN 

IF  ( UIN  1 UI  .NE.  1.0)  BZ=PBZ  ( UIN  1 UI  , P 1 U,  NU  WBEN) 

ZI=  (FAZ  (UIH1UI.P1B,  NUWBN)  *ZIH  1 ♦ BZ  *F2B  *DELX*  ( 1 . +RR  AT  1 PN)  /2  . ) 

$/R  SAT  1 P N 

ZB=  (ZI*ZIN1)/2.  S D2=PD2(N,ZI,C1) 

...EVALUATE  FUNCTIONS  TO  OBTAIN  A NEk  VALUE  FOB  HSTA8  AT 
...STATION  I 

RD2B  = FRD2  (N.ZB.C1B)  $ D2B=  F E2  (N  , Z B,  C 1 B) 

IP (FLOW  . EQ.  "TURBULENT")  GO  TO  340 
CDB=  FC0L (HB,B,G, NDB, THETB, W,RC2B,D2D2UB) 

GO  TC  360 

340  IP  (NCOS  . EQ.  1)  GO  TO  350 
D2UB=D2B/D2D2UB 

A13=FA1  (H12B)  $ RD2UB=B  D2B/C2  02  UB 

CFIB=2. *A1 B/RD2UB*»N  $ CPB= FCP ( AB, BD2B , D2D2UB) 

J8=FJ  (H12B.CFIB)  $ JSPB  = P JSP  ( A 1b , H 1 2B,  HB) 

PIB=FPI  (H12B,  D2UB,  CFIB,  DNSDXNS) 
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JP B=  PJP  (JB.PIB,  KCDS) 

IP  (KCDS  . EQ.  1)  GO  TO  350 

CDB*  PCDTP  (CFB,  RD2US.N,  JPB.JB,  JS  PB,  MB.  H 1 2 B,  PIB ) S GO  TO  360 

350  CDB=PCDTRT(D2D2UB,RD2B) 

360  F 4 B=  P F 4 (BD2B,M,CCB,AB,  D2D2UB,  USB) 

IP (ERROR  .AND.  PLOW  . EQ.  "TURBULENT")  WRITE(6,923)  A1B,CPIB,JB, 
JJSPB.PIB, JPB,CDB,F4B 
BH= 1.0 

IP(UIM1UI  .NE.  1.0)  BH=PBd(UIM1UI,P3B) 

HIS*  PAH  (UIM10I,  F3B)  *HIH  IS*  BU*  F4  B/ZB*  DKLX 

...ALLO*  POR  INTERVAL  SUBDIVISION  IF  NECESSARY 

I P ( ICOU NT  .GT.  1 .OR.  I .GT.  1)  GO  TO  370 
LI  M2=LFIX  (ABS(UIS-HlM1S)/0.  00  25)  *1 
LIM=HAX0(LIM1,LIH2)  $ LIMUP=50 
LI M=  MI  NO  (LIM.LINUP) 

IF  (P  LON  .EQ.  "LAMINAR"  .AND.  KPOINT  .LE.  10)  LIH=20 
IF  (KPOI NT  .LE.  2 .OR.  FLON  .NE.  PLOWIN 1 ) LIH=50 
IF  (LIM  .GT.  1)  GO  TO  230 

...EVALUATE  FUNCTIONS  TO  FIND  A NEU  VALUE  OP  CUT  FOR  HTCOR  R=.  T SUE. 

370  IP  (.  NOT.  (HTCORR) ) GO  TO  390 

HSP= (HIS-HIM1S)/DELX  $ DD  2DXD2*  (D2-D 21 M 1)  / (D2B* DEL X) 

P58=PF5 (DMSDXMS.HSB, D1D2BfG,MDB,HSP,DD2DXD2) 

P6B=PP6 (HSB.BP, BB, DHSDX MS , D 1D2B , G , MDB) 

I?  (F5B  .NE.  3.3)  GO  TO  380 
CHT=CHTIM1*PbB*DELX  $ GO  TO  390 
380  r.'!T^F6B/F5B*  (CHTIM  1-F6B/F5B)  *BXP(~P5B*DELX) 

...A  NEW  VALUE  FOR  HSTAR  AT  STATION  I IS  NOW  KNOWN.  ITERATE 
...TO  FIND  THE  CORRESPONDING  H VALUE 

39)  11=3  i HI=HINNER(1)=PH(HIS,PSIB) 

IF  (ERROR)  WRITE  (6, 924)  H INNER  (1) 

4)3  LI  = I It  1 

IF  (H I .LE.  1.515)  GO  TO  510 

I F (H I .LE.  2.3)  GO  TO  413  $ W R ITE  (6 , 922 ) X I $ GO  TO  633 

4 10  I?  (FLOW  . EQ.  "TURBULENT")  GO  TO  420 

D1UD=FD1UDL  (HI)  t G1=PG1L(HI)  * GO  TO  430 
420  D1  UD=FD1  UDT  (HI)  $ G1=FG1T(Hl) 

433  IF  (ABS(TUETI)  . LE.  2.3)  GO  TO  440 
PSI=1.0  I GO  TO  450 

443  PS  I 3P=PPSI3P  (HI,  THETI)  $ PS1 1 2*  FPSI 1 ? (D  1 UD,  HI  ,TH  ETI  ,G  1 ) 
PSI=FPSI  (PSI12,MDI,PSI0P) 

45  ) HINEW=H INNER  (II*  1)  =PH  (H1S,PSI) 

IP(ERROR)  W RITE  (6,925)  II  ♦ 1 , BIN  NEB  (II  ♦ 1 ) 

IP  (ABS  ( (H1-H1NEW) /HINEW)  .LE.  EPS/10.)  GO  TO  470 
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PBOGRAH  CON  PBL  . 


IP  (II  .GE.  20)  GO  TO  460 
KISHINEV  $ GO  TO  400 

460  UBITE  (6,926)  XI , (I  P 1 , HINN ER  (IP  1 ) , IP  1*  1, 2 1 ) S GO  TO  600 

470  HI=HHAIN  (I*  1 ) = (Sl  + HlNEW)/2. 

IP(EBBOB)  WBITE(6,927)  I*  1 , UN  AI  N (I*  1 ) , BI , CHT 
C 

C...THE  VALUE  OP  d AT  STATION  I POB  TB1S  ITEBATION  IN  THE  NAIN 
C.  • • LOOP  IS  NON  KNOWN.  CONTINUE  ITERATING  IN  NAIN  LOOP  IP 
C...NECESSARX 

C 

IF  (ABS  ( (HI-HHAIN  (I)  )/Hl)  .LB.  EPS)  GO  TO  500 
IP  (I  .GE.  20)  GO  TO  490 

IP  (.NOT.  (ICOUNT  .EQ.  1 .AND.  I .GE.  5 .AND.  KTBIP  . EQ.  0)) 

I GO  TO  480 

KT  RI F=  1 $ LIH=  20  $ GO  TO  230 

480  HIOLD=Hi  $ GO  TO  2B3 

490  WRIT E (6 , 92  8)  XI , (1P2,  HH  AI  N (IP2)  ,1  P2=  1 ,2  1)  S GOTO  600 

C 

C... CONVERGENCE  FOR  H AT  STATION  I IN  THE  NAIN  ITEBATION  LOOP  HAS 
C...OCCUBRED.  CHECK  FOB  SEPARATION 
C 

500  HI= (HI+HNAIN  (I) )/2. 

IP  (HI  .GT.  1.515)  GO  TO  520 
510  IP  (FLOW  .EQ.  "LANINAB")  WRITE  (6,929)  XI 

IP  (FLOW  .EQ.  "TURBULENT")  WRITE  (6,930)  XI 
GO  TO  600 
C 

C... CALCULATE  AND  PRINT  OUTPUT  VARIABLES  OP  INTEREST 

C 

520  IP  (HI  .LE.  2.0)  GO  TO  530  S WRITE  (6 ,922)  XI  $ GO  TO  600 
533  IP  (FLOW  .EQ.  "TURBULENT")  GO  TO  540 

A=FAI(HI)  $ H 12=PH12L  (HI)  $ D 1UD=PD1  UDL  (HI) 

GO  TO  550 

540  A=  PAT  (HI)  S H 1 2=PH  1 2T  (HI)  t D1UD=PD1UDT(HI) 

550  D2D2  0=PD2D2U  (R,  G,  N DI  ,HIS,  THETI) 

D1D2=FD1D2  (H  1 2 , D2D2U,  R,  G,  NDI,  HIS,  THETI,  CHT) 

560  D2=FC2  (N,ZI,C1)  $ R D2«PBD  2 ( N,  ZI , C 1) 

RD2U=FBD2U  (R D2 , D2D2U,  R , G,  NDI,  THETI)  i D2U=D2/D2D2U 
D1=D2»D1D2  $ CF=FCP  (A  ,RD2  ,N  .D2D20)  $ RCORR=RI*D1 

QDIN=FQDIH(CF,S,BI,CUT,G,NDI)  S D999= (H 1 2 *D2) / (D1UD*D2D2U ) 
LCCUNT=LCOUNT+1 

IP  (LCOUNT  .LT.  LIN  .AND.  (.NOT.  ERROR) ) GO  TO  580 
IP ( ERROR  .OR.  (IP  .NE.  17  .AND.  FLOAT(IP-17)  . NE.  25.*PLOAT({ 
$IP-17)/25) ) ) GO  TO  570 
WRITE  (6,901)  S WRITE  (6,  9 1 4) 

WRITE(4,901)  t WRITE  (4, 91 5) 

570  IP=IP+1 

WRITE  (6,931)  XI,ZI,HI,D1D2,D1,U2,D999,RD2,CP,QDIH 
WRITE  (4,919)  XI  , BCORR 
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P BOG  BAM  CORPBL  . 


...CHECK  FOR  TRANSITION 
580  FLORIN Is  PLOW 

IP (P LOU  .EQ.  "TURBULENT"  .08.  NOTBAN  .OR.  (XTBAN  . NE  . 0.0  .AND.  II 
S . II.  XTBAN) ) GO  TO  59J 
PPARI=FPPAR  (C1,DMSDXMS,D2U,  R,  G,  HDI , TH  ETI) 

IF  (ZSTART  .NE.  0.0  -AND.  KPOINT  . EQ.  2)  PP ABI M 1=FPP AR (C 1IH 1 , 
IDMSOXMS, D2UIM1,H.G,MDIN1,TliETIM1) 

IF  (XTBAN  .NE.  0.0  .AND.  XI  .GE.  XTRAN)  KTRA  NS=  1 
CALL  IRA  NS  (Z  1 , HI) 

...SHIPT  I QUANTITIES  TO  1-1  QUANTITIES  AND  READ  THE  NEXT  POINT 
...FOR  THE  BOUNDABX  LAYER  COMPUTATIONS 

590  DHDXIM1=0.0 

IF  (KPOINT  .GT.  1)  DHDXIH1=(HI~H1N1)/DBLX 

XIM1  = X1  $ MDI M 1=MDI  $ MLIM1S=MDIS  $ DELXIN  1=DELX 

TO  ETIM  I^THETI  $ TNTIIH1 =T WT II  * RIM1=RI  $ YIN1=YI 

H1.11  = HI  $ HIM  1S-HIS  $ Z1N1=ZI  S PPARIM1=PPARI 

BI M 1-B1  $ D21M  1=D2  $ CHTIN1=U1T  $ D2UIH1=D2U 

C 1 1 M 1 = C 1 

IP (LCOUNT  .LT.  L1M ) GO  TO  240 
GO  TC  210 
600  STOP 


FORMAT  STATEMENTS 


9)0 

901 

902 

903 

904 

905 

906 

907 

908 

909 

910 

911 

912 


FORMAT  (8A10,/,4A13) 

FORMAT ( 1 H 1 ) 

FORMAT (/,43X,*INPUT  ERROR: 
i *" TURBU  LFNT " * ) 

FORMAT  (/,40X, “INPUT  ERROR: 
i*"AXISYM"») 

FORMAT  (/, 4 )X,  “INPUT  ERROR: 
$*"MSTAB"») 

FORMAT  (/,40X,*INPUT  ERROR: 


FLOW  MUST  EQUAL  "LAMINAR"  OB  », 
GEOM  MUST  EQUAL  "PLANE  2-D"  OR  «, 
MTYPE  MUST  EQUAL  "MACH"  OR  *, 
TTYPE  MUST  EQUAL  "THETA"  OR  *, 


$»"TNTINF"») 

FORMAT  (/, 4 )X,  "INPUT  ERROR:  iiSTART  MUST  BE  .GT.  1.515  AND  .LE.  ", 

i"  2.0") 

PORMAT (/,40X, "INPUT  ERROR:  BSTART  MUST  BE  .GE.  -0.1988  AND  ", 
$".LE.  2.0") 

PORMAT  (/,J0X, "INPUT  ERROR:  ZSTART,FSTINT, N MUST  BE  .GE.  0.0  ", 

$ "A N D M1NF,RE1NFL,RL,>T,SL,ST,EPS  MUST  BE  . GT.  0.0") 

FORMAT (/,40X, "INPUT  ERROR:  G MUST  BE  .GE.  1.0  AND  . LE.  1.67") 
FORMAT  (/,40X, "INPUT,  ERROR:  NINF»  MUST  BE  .LE.  SQRT  ( (G“  1 . ) /"  , 

*" (G- 1 . ) ) " ) 

FORMAT  (//,5dX, "EXECUTION  DELETED") 

PORMAT (32X, "COMPRESSI BLE  BOUNDARY  LAYER  RESULTS — COMPUTATIONAL", 
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PROGRAM  CCHPBL  . 


**•  METHOD  II  OF  H A LZ",///, 26X, 8A  10 ,/ .58X , "INPUT  PARAMETERS:",//, 

$36  X,  "FLOW  =". 1X,A10,6X."GEOH  ="  , IX  , A 1 3,  5X  , "NT  YPE  =",1X,A10,/, 

$35  X,  "TT  Y PE  «". 1X,A10,4X,"ZSTABT  = , G 1 0.  3, 5X  , "HSTART  =",G10.4,/, 

$ 34  X#  "BSTART  ="  ,G  1 3. 4,  7X.  " HI  Nr  = ",  G 1 3.  4,  SX,  "REI  NFL  =",G10.4,/, 

$J5X, "XTBAN  =",G10. 4,5X,"FSTINT  =»,G 1 3. 4, 10X, "G  =".G10.4,/, 

$ 3 9 X , "W  = ",GU.4.9X,"RL  =",G10.4,9X,  "RT  = ",G13.4,/, 

$ 38  X , "S  L — " ,G10.4,9X,"  ST  = ",  G1  3.  4,  8X,  " EPS  =",G10.4,/, 

$ 34  X,  "NOT  RAN  = " , IX  , L 1 , 1 3X,  "HTCORB  =",  IX , L 1 , 1 4X,  "ERROR  =",1X,L1, 
$///, 62X, "RESULTS: " ,/) 

913  FORMAT  (//,  5 3 X , "CORR  ECTED  BOUNDARY  CO-OR  DI  NATES"  ,////,  2 5 X,  8A  1 0 ,///) 

914  FORMAT  (4 X,  "AXIAL",  1 OX,  "Z  = ",  11  X.  "H^",9X,  "SHAPE",  5X,  "DISPLACEMENT", 
iJX, "MOMENTUM", 4X, "99. 9*  B. L . " , 2 X, "MOM . THICK. ", 3X ," LOCAL  SKIN", 
P2X, "DIMENSIONLESS",/, 3X , "LOCATION", 3 X , "D2« (RD2**N) ",4X, 

i" (D3/D2) U", 6 X," FACTOR", 5X, "THICKNESS", 4X , "THICKNESS ",  4 X, 

$" THICKNESS ",5X, "PEYNOLDS" , 5X, "FRICTION", 4X, "LOCAL  «ALL",/,5X, 

$"[  L ]" , 1 0X , "[  L ]"  , 2 1 X , " (D1/D2)  " , 8X,  "[  L ] ",  1 OX , "(  L 1 0X,  "[  L ]",  9X  , 

$ "NUMBER  ",  4 X,  "COE'FFICI  ENT",  4X,  "HEAT  FLUX",//) 

915  FORMAT (34X, "AXIAL  LOC AT  ION" , 29 X , "AXIS YM METRIC  RADIUS,  OH",/, 
$«0X,"£L]".31X,"2-D  DISTANCE  PROM  CENT  BH  LINE",/,  87X,  »[  L ]"  ,//) 

916  FORMAT  (5F10.  5) 

917  FORMAT  (3  ( IX,  G 12.6)  ,5X,  A4,4X,4  (1X.G12.6)  ,2  (5X,  A4,4X)  ,/) 

918  FORMAT  (3  (1X.G12.6)  ,4(5X,A«,4X)  , IX ,G  1 2.6 , 2 (5X,  A4,4X)  ,/) 

919  FO  RM  AT  (2  ( 35X  , G 1 2. 6)  ,/) 

920  FORMAT  ( 17X,  "11HAIH  ( 1)=",G12.6) 

921  FORMAT  (22X,"HBAR=",G12.  6) 

922  FORMAT<1X.G12.6,//,24X,”PHOBABLE  CONVERGENCE  PROBLEMS:  H ", 

$ " EXC  ELDS  2.0  — BOUNDARY  LAYBR  COMPUTATIONS  TERMINATED") 

923  FORMAT  <2X , "A  = ",G  12.  6,  1X,"CFI=", G12.  6,  IX,  "J=",  G12.  6,  IX, 
$"J*'=»,G12.6.  IX,  "P1=",G12.6,  IX,  "J'  = ",G12.6,  1 X,  "CD="  ,G  1 2. 6 , 

$ 1X,"F4=",G1Z.  6) 

924  FORMAT  (16X,"HINNER ( 1)=",G12.6) 

925  PORMA1  ( 1 bX,  "HINNEB  (",  12  ,")  =" , G 1 2. 6) 

926  FORMAT (IX. G12. 6,//, 38X, "CONVERGENCE  PROBLEMS  FOR  H GIVEN  HSTAB: ", 
*"  EXECUTION  HALTED",//,  (54X  , "HI  NN  ER  ("  , 12  . ")  =",  G 1 2 . 6) ) 

927  PORMAT  ( 1 7X,  "HNAIN  (",I2,  ")  =", G 1 2. 6 , 2 OX  . "B=", G1 2.  6,  20X,  "CHT=  ",G  12.  6 ) 

928  FORMAT (IX, G12. 6,//, 32X, "CONVERGENCE  PROBLEMS  FOR  H IN  MAIN  ", 

$ "I  TERATICN  LOOP:  EXECUTION  HALTED »,//,( 55 X,  "HH AIN  (" ,12  ,")  =", 

$G1  2.  6)  ) 

929  FORMAT(1X,G12. 6, //,37X, "LAMINAR  SEPAR ATION— BOUNDARY  LAYER  ", 

$ "COMPUTATIONS  TERMINATED") 

933  FORMAT ( IX. G12. 6, //,38X, "TURBULENT  S EPAR ATION  — BOU  NDAR Y LAYER  ", 

$ "COMPUTATIONS  TERMINATED") 

931  FORMAT (10 (1X.G12. 6) ,/) 

END 
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FUNCTION  PJP  . 


FUNCTION  FJP  (J.P.KCDS) 


C. 

C. 

c. 

c 


FUNCTION  SUBPROGRAM  FJP  EVALUATES  J'  FOR  FELSCH'S  TUKBULi-NT 
CISSIPA1ION  LAW  IF  J AND  P ARE  IN  THE  PROPtB  RANGES.  OTHER- 
WISE KCDS  IS  SET  TO  1 AND  THE  HCTTa-T RUC KFNBRODT  DISSIPATION 
LAW  IS  USED  IN  TH E MAIN  PnOGBAM 

IMPLICIT  RSAL(J) 

DINENSION  JP  (17) 

DATA  JP/-300. ,-200. ,-100. ,-50. ,-40. ,-20. ,0. ,20. ,40. ,60. ,83. , 

tin. ,-jjj.,-20j. ,-i jj.,  ). ,ioj./ 

JP  10  0 (X)  = -9.  j 188658 708729+1  . ISO 34 54  1 3206 3«X 
S + J.55  5yJ15  0Jt2227t'-J*X*  + 2-J.  183  43  543J)5468ii-3*X**3 
$+0.1«1t0S1333171fcL-5*X**4 
JP80  IX)  =-4.  34  1373  17o1d76+  1.  30  6B42b01 158  4*X 
S-0.  1 12  500671  59924E-  1»X*  *2  *-0.  472 755 1 0 54b  JOBE-  4 ♦ X **3 
J?6 J (X)  = - 1 . J 3 3 88478  48 44 3 + 1 . 35o 290 73 1 ) 73 5* X 
S-0. 149 025228 136 89E-1*X* *2  +0. 76509 Ob 55 76  344 E-4*X** 3 
JP4J(X)=1.5432J4511719)+1.b546l35905b38*X 
$-0.34 08 3 01b059t81£-1*X* *2+0. 443 73331872342 E-3*X**3 
,p-).223769175a7)56E-S*X»*4 
JP20 tX>  =3. 607 5 10531 179 4*1.8 5352 83201 0 70  + X 
X-  3.  4 7 11  u 1955  26  /22t- 1 ♦ X * » 2 ♦ 0 . 7 0 )32  7 3 1 7 58 5 16  E-3*  X**3 
S-0.384v43B69G759bi.-5*X**4 
JP  3 (X)  =7.  ) 1824 3b 39602b*  1.  84 06b  1 19  17541  *X 
X-O.b It  1 12228a 91 89E-1*X» *2+0.81749 4 JOB 389y 9 E-3*X**3 
S-  3.4  68  43)30978392  E- 5*  X**4 
JPM20  (X)  = 6.  9 2 3b  17  2 182  76 8*  2.  04  55 3 196  25 133* X 
S- J.  o55j6  35  7 7 1 82  )'»E-  1*X**2+  3.  1 03  37  1753  1S598E-2*  X**3 
S-0. 6 40 2 58868  0 5832 L- 5* X*«4 
JPM4  3 (X)  = 1 ).  75d89  36  384  9 9 + 2.  )b  23  62  734  6 909»  X 
S-0. 6 9J807  021 00035E-1*X**2  + 0. 1 1 7SC64  11  1557  1b-2*X**3 
S - 3 . 6 99  3 25  1 2d  1 d 3 1 2 --  5*  X * * 4 
JPN5  0(X)=  11.  7 34292 52892  1*2. 0573 0884b 7056*X 
S-J.b  9739322792  J3+i.-  1*  <**2  + 3.  1 19  353292  19484£-2*X**3 
S-0.  7 089C4  356  17  54  36-5*1*  *4 
JPM1 30 jX)  =1  3.  2921  3750 33 2) *2.  1 79 898 55  Ibdb 1 * X 
S-). 751241334006462-1* X* *2 +0.126327418746 06 E-2*X**3 
S-).  74244  )9  5 J 9 41bO_-5*X**4 
J P 8 200  (X) - 14. 40849 2394423  + 2. 4 8 3 9430984 792*X 
J. 8 6 1 J 38  58 7833396 -1*X**2+J.  15  35 34  3 3958 791 S-P*X**3 
S-0. 396720735341 35E-5*X**4 
J?.M33)(X)  = 15. 45615336  9 6 33 +2. 59849117  o 7399*X 
.+  -3.  3 702030961  1832.-,- 1*X*»  2 + 0.  14480538b  16  36  3S-2*X**  3 
S-  >. bUy  155243651 8 b..- 5*  X**4 
JP  1 3 1 ( a)  =47.61  + 3 . 4 * X 
JP  J 1 (X)  - -i3.9)+J.4*X 
JPN  1 01  ja)  =52.82  + 0.  4 * X 
JP3201  (X)  =57.40+  3.  4«X 
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FUNCTION  rJl  . 


JPfIJOl  (X)  =6J.79*0.4»X 

PJP 

500 

JUNO 

= JPN300  ( P)  i JLINL  = JP1  00  (?) 

FJP 

510 

1 F (P 

• GT.  70.)  JLIMU=0?NJ01  (P) 

FJP 

520 

1 3 <P 

-GT.  70.)  JLINL=J?101  (P) 

FJP 

530 

1P(.NGI.  (P  .LT.  0.0  .OR.  J . OT.  OHIO  .OR.  J .LT.  JLINL) ) 

GO  TO 

10 

FJP 

540 

KCDS  = 

1 S RETURN 

FJP 

550 

1 ) 

if  { r 

. UT.  7 ).  0)  GO  TO  1 2) 

FJP 

56  0 

IF  (J 

. LE.  JPH200  (P)  ) GO  TO  20 

FJP 

570 

1=  1 

•3  FRA  C=  (J  - 0 PN  )J)(P)  )/(.iri2  ) 3 ( P)  - J PN  3 0 3 (P)  ) 

$ 

JO  TC 

16;  0 

FJP 

58  0 

20 

IF  (J 

• L*-.  JPM100  (P)  ) GO  TO  30 

FJP 

59  0 

1-2 

J.  EVAC*  (.J-JPN2  ) )(?)  ) / (JIM  1 ))  (?)  -JPN200  (P)  ) 

A 

GO  TC 

1b  0 

FJP 

6)0 

30 

IF  (J 

. Lt.  JPMSJ (P) ) GO  10  40 

FJP 

6 10 

1=  J 

i i'RAC=  (J-JPN1  ) )(?)  )/{JIK5)  (?)  -JPS1  0)  (P)  ) 

£ 

GO  TO  160 

FJP 

62  0 

40 

IF  (J 

. LE.  JPH40  (P)  ) GO  TO  50 

FJP 

630 

1 = 4 

i FRAC=  (.W?N5)  (P)  ) /(J?!140  (P)  -JPM5)  (P)  ) $ 

CiO 

TO  160 

FJP 

64  0 

50 

IF  (J 

.L_.  J ? F 2 0 ( I ) ) GO  TO  60 

FJP 

650 

1 = 5 

4>  FRAv:=  (G-JPN40  (P) ) /(JP»,2  )(?)  -JPN40  (P)  ) i 

GO 

TO  16  0 

FJP 

66  0 

6u 

1 F ( J 

. LF.  ,TPC  (P)  ) GO  TO  70 

FJP 

670 

1=  6 

* PKAC'=  (0-JPN2J  (?)  ) / (JPJ(P)  - JPH20  (P)  ) $ 

GO 

TO 

16  0 

FJP 

6U  0 

70 

IF  <J 

-LE.  JP20  (?) ) GO  TO  BO 

FJP 

690 

1=7 

i FKAC=  (J-J?3(P))/(JPx)  (P)-J?0(P)  ) $ GO 

TO 

160 

FJP 

700 

BO 

Ir  (.! 

• Li.  .1 P 4 0 ( P ) ) GO  TO  90 

FJP 

7 10 

1=  .) 

!>  FwAC=  (J-JP2  )(?))/ (JP4  ) (P)  -J?2  1 (P)  ) $ 

GO 

TO 

16  0 

FJP 

720 

9 3 

IF  (J 

.LG.  JPtO(P))  GO  TO  100 

FJP 

730 

1 = 0 

a F«AC=  (0-JP4  )(?))/ (JP6  0(P) -JP4J  (P)  ) $ 

GO 

TO 

16  ) 

FJP 

74  0 

IT) 

IP  (J 

.LF.  JitiO(P)  ) GO  TO  110 

FJP 

750 

1=  10 

F1AC=  (J-J?r.  ) (P)  ) /(JPdJ  (?)  -JP60  (?)  ) i 

GO 

10 

16  0 

FJP 

760 

in 

1=  1 1 

a FhAC=  (J-JPoO  (P))/ (J?  100  IP) -0P80  (?)  ) £ 

GO  TO  160 

FJP 

770 

12) 

I F ( J 

. LF.  J PN2)  1 (?)  ) GO  TO  100 

FJP 

780 

1=  1 ) 

a Fk.i^=  (J-JPMJOI  (P) ) / (JPH20  1 (P)  - J PM30 1 (P)  ) 

FJP 

790 

go  rc 

i t>  ) 

FJP 

BOO 

1 ) ) 

IF  (.1 

. La.  JPNIOI(PJ)  GO  TO  140 

FJP 

810 

1=14 

i Fr'AC=  (J-JPK2  01  (P)  ) /(JPN1  0 1 (p)  -JPN20 1 (P)  ) 

FJP 

82  0 

GO  TC 

1 60 

FJP 

83  0 

140 

IF  (J 

.Lt.  Jr 0 1 (P)  ) GO  TO  15  0 

pjr 

840 

r=  15 

i ?RAC=  (J-JPN  101  (P)  ) / (J?01  (P)  -JP.M101  (P(  ) 

5 

JO  TO  160 

FJP 

850 

150 

1=  16 

v FT  AC=  (J-JP01  (?) ) /(JP13  1 (P)  -JPO  1 (P)  ) 

FJP 

86) 

16) 

FJ?=JP(r)  ♦ FiiAC*  tJP(I*  1 ) -JP(I)  ) 

PJP 

070 

hF.ru  r. 

i; 

FJP 

880 

L'fD 

FJP 

890 

77 
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SUBROUTINE  TRANS  . 


SUBROUTINE  TRANS  (L  , ti ) 

C 

C. . . SUBROUTINE  TRANS  CHECKS  FOR  TRANSITION  FROM  A LAMINAR  TO  A TURBU- 
C...T.ENT  BOUNDARY  LAYER  USING  THE  REYNOLDS  NUMBER  BASED  ON  THE  INCCM- 
C.  . . PRISSIBLE,  ADIABATIC  MOMENTUM  TfiiCK  NESS.  THE  POINT  OF  INSTABILITY 
C...IS  FOUND  USING  A CURVE  FIT  TO  THE  INCOMPRESSIBLE  COMPUTATIONS  OF 
C...WAZZAN,  ET.  AL.  THE  TRANSITION  POINT  IS  1HEN  FOUND  USING  A COR- 
C...  RELATION  TO  liRANVUl-'S  DATA  FOR  THE  DISTANCE  BETWEEN  THF 
C. .. INSTABILITY . AND  TRANSITION  POINTS,  AND  A CORRECTION  FOR  FRRESTKEAM 
...TURBULENCE  EFFECTS. 

IMPLICIT  RFAi.  (J,  K,  N) 

REAL  INTGPAL 

CO MMCN/CTRA NB/RD2 , N,  R , R I,  S,  ST,  FLO  W ,G  , M TYE  E ,TT  YPE, MINF, KTRANS, 
iK  INST,  RD2U,  PI  A R I,  P PAR  I Ml,  DELX,  F ST  I NT 
IF  (KTRhNS  .DQ.  1)  SC  TO  23 
IF  (KINSI  . E<*.  1 ) GO  TO  1 J 

...CHECK  FOR  INSTABILITY  USING  WAZZAN’S  RESULTS 

h02INSI=10.**  (-5077986.  1 033 98  2*  15932  1 0 1 . 403  096*ii 
$- 199842ob. 72,.  1 77* H* *2*12526914. 450584*8**3 
I- 3924)56. 68  349  J9*H»*  4*491  4 15.  2573  57  8 5*11**5) 

IF  (H  .81.  1.  625)  l!D2iNST  = 5093.  *177209.  * (H-1. 625) 

IF  (PD2U  .LT.  KD2i.NST)  RETURN 
WRITE  (6,990) 

IF  (FSTi  NT  .GT.  2.  16  3)  GO  TO  20 
KINST-1  i DEI.SUM=0.  0 $ INTGRAL=0.0 

TCiffF-  (900.-760.  * FSTI  NT  * 1 59  . *FSTI  NT*  * 2)  /900.  $ RETURN 

C 

C...FIWD  THE  TRANSITION  POINT  USING  A CORRELATION  TO  GRANVILLE'S 
C . . . DATA 
C 

1 ) I NTGRAL- 1NTG RA L* ( P PARI* PP AR IM 1) /2. *QLLX 

DcLSU.M-  DELSUJi*DELX 
FPAPM=INI3RAL/DELSUM 

HD2T  RA  N-  R D2  IN  SI*  (450.  *4  00.*  EXP  (60 . *PPARM)  ) *TCORP 
I f (H  C2U  .LT.  D02T  RAN)  RETURN 
2 0 WRITE  (6,901) 

N=).26o  S h = hV  A S=ST  t E'  LOW  = "TURBULENT" 

Z = Z*  (c.  02* * (-0.7  12)  ) 

UP.  10  ;>  i 
3 

C FORMA  I jTATEMENTG 

c 

900  FORMAT  (/  , 5 3 X , " L \ M I NA  R INSTABILlTY  POINT",//) 

901  format  (/,r.2x, "lami nap- turbulent  transition",//) 
fn  r 


TRA 

10 

TRA 

20 

TRA 

30 

TRA 

4) 

TRA 

50 

TRA 

60 

TRA 

70 

TRA 

BO 

TRA 

90 

TEA 

100 

TRA 

110 

TRA 

120 

TRA 

130 

TRA 

140 

TRA 

150 

TRA 

160 

TRA 

170 

TRA 

18  ) 

TRA 

190 

TRA 

2)0 

TRA 

210 

TRA 

220 

TEA 

230 

TPA 

240 

TRA 

250 

TPA 

260 

TRA 

27  3 

TRA 

280 

TRA 

29) 

TRA 

300 

TPA 

310 

TRA 

320 

TRA 

330 

TRA 

340 

TEA 

350 

TEA 

36  0 

TRA 

370 

TRA 

38  3 

TRA 

390 

TRA 

40  ) 

TRA 

413 

TRA 

42  0 

TRA 

430 

TPA 

44  0 

TRA 

450 

TRA 

46  0 

TRA 

473 

TRA 

483 

78 


non  r-.  r n 


p 


THIS  PAGE  IS  BEST  QUALITY  PKACJICABIig 
FROM  COPY  FURNISHED  TO  DOG  r - 


APPENDIX  A.  COMPuESSaULE  BOUNDARY  LAYER  PROG  RAM. .. COMPBL  PAGE  A-16 

subroutine  nsort  . 


SUBROUTINE  NSORT  IM, MS) 

HSO 

10 

c 

HSO 

20 

c. . . 

SU3ROUI  INF.  NSORT  DETERMINES  HHETHER  M IS  A MACH  NUMBER  OR 

HSO 

JO 

c... 

NSTA  H AND  CONVERTS  fl  TO  A MACIJ  NUMBER  AND  MS  TO  AN  MSTAN 

HSO 

40 

c 

HSO 

50 

IMPLICIT  RSAL(J,M,N) 

HSO 

60 

CCMMON/CT  RANS/RD2,  N,  R,  NT,  S,  ST  , PLO«  , u , IT  Y PE,  TTY  P E,  MINE,  KTRANS, 

HSO 

70 

GKINST,hD2U,PPARI,PPAKIM1,DKLX,PSTINT 

HSO 

80 

c 

HSO 

90 

c.. . 

STATEMENT  FUNCTIONS 

MSO 

100 

c 

HSO 

110 

EM  (G  ,MS)  =SioT  (2./  i <S ♦ 1 .)  »NS**2/  ( 1.-  iG-  1.  )/  (G+1.  ) *».S**2)  ) 

MSO 

120 

FIG  (G,  M)  =SUR'i  ( (G*  1.J/2.  *’1**2/  (1  .♦  (G-1  .)/2.*H**2) ) 

MSO 

1 JO 

c 

MSO 

140 

c. . . 

DO  SORTING  AN J CONVERTING 

MSO 

150 

c 

MSO 

160 

IF  IMTYP  E . EU.  "AST  Ah")  uU  TO  10 

MSO 

170 

MG  - F MS  (G  , 1)  A RETURN 

MSO 

180 

10 

MS=M  J>  M = FK  (G  , M) 

MSO 

190 

K’G  TURN 

MSO 

200 

END 

HSO 

210 

t 


onr  non 
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SUBROUTINE  TSORT  . 


SUBROUTINE  TSORT  (I  i.ETA  , TWTIHP  , 3 . ND) 

C 

C.. .SUBPOUTINL  TSORT  DETERMINES  WHETHER  THE  INPUT  TEMPERATURE  DATA 
C...IS  IN  TERMS  OF  THE1A  OR  TWALL/TINF  AND  CALCULATES  AND  STORES 
C... THETA  VALUES  . N THETA,.  TWALL/TINF  yALUES  IN  TWTINP,  AND  b VALUES 
C..  .IN  D 
C 

IMPLICIT  RFAL(.J,N,N) 

C0'1MCN/CTRANS/RB2,N,K,KT,S,ST,FL0K,G,MTYPE,TTYPE,M1NF/KTHANS, 
$K1NST,RD21J,PPARI,  PPARIM  1 , DELY  , P STINT 

...STATEMENT  FUNCTIONS 

FTHKTA  (R.G.MD,  MINF ,TWTI NF) = ( (1.  *n*  (G-  1.  ) /2.  *ND**2)  - (1 . ♦ (G-  1 . ) /2. 
f*ND**2)  / (1.  ♦ (G-1  .)  /2.*MINF**2)  *TWTI  NP)  / (R*  (G-1.)/2.*MD**2) 
FTWT1NF  (R,G,  KD,  Mi  NF,  THETA)  = ( 1 . » (G  - 1 .)  /2.  »MI  NF*  * 2)  / ( 1 . ♦ (G-  1.)  /2. 
i*MD**2)  « (1.*R«  (G-1.)/2.  *MD**2*  ( 1.  -THETA)  ) 

FB  (THETA  , R,  o,  MD)  =TH ETA*  U*  (G- 1 . ) /2  . *MD**  2 

...DO  SORTING  AND  CONVERTING 

IF ( T T Y P L -EG.  "TWTINF")  GO  TO  10 

TW II NP  - FT  WT  i N F (R  , G , KD  , M IN  P,  TH  IT  A)  $ GO  TO  20 

10  TWTINFsTULTA  S THETA=FTHETA  (R,G,M  J,  MINF,TrfETA> 

2 0 ['=  FR  (Tl.fc  IA  , R , G,  ME) 

RETURN 

END 


TSO  10 
TSO  20 
TSO  3 0 
TSO  40 
TSO  S3 
TSO  60 
TSO  7 3 
TSO  80 
TSO  90 
TSO  100 
TSO  110 
TSO  120 
TSO  130 
TSO  140 
TSO  150 
TSO  16  "> 
TSO  17 7 
TSO  180 
TSO  190 
TSO  200 
TSO  210 
TSO  220 
TSO  2 J 0 
TSO  240 
TSO  250 
TSO  260 
TSO  273 
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SUBROUTINE  IE8ROH  . 


SU  ERCUTIN  E IRRROR  (DELX  , MDI , (1DIS  ,G  , H I , GEGH  , TWTI N P)  .RETURNS  (I) 

C 

C. .. SUBROUTINE  1ERROR  CHECKS  FOR  INPUT  ERRORS  ON  THE  VALUES  OF  XI, YI, 
C.. . HDI, HDIS.RI,  AND  TNTINP 

C 

IMPLICIT  R2AL(J,H,N) 

K=0 


IF 

(DFLX 

. GT . 0.0)  GO 

TO  10  i 

K=K*1  S 

WRITE  (6, 

900) 

1 ) 

IF 

(MDI  . 

GT.  0.  3)  GO 

TO  20  i 

K=K*1  $ 

WRITE  (6 

,9 

01) 

20 

IF 

(MD1S 

• L E.  S v RT  ( (G 

♦ 1.)/(G-1.)  ) 

•AND.  MDIS 

• GT.  0. 

0) 

bO 

TO 

K = 

«♦  1 

S WRITE  (6, 

9 02) 

30 

IF 

(.NOT. 

(GECM  .Ew.  " 

AXISYM"  .AND 

. HI  .LE.  0. 

,0))  GO 

TO 

40 

K = 

K*  1 

i WRITE (6, 

90  J) 

40 

IF 

(.NOT. 

(GEOM  .Ew  " 

PLANE'  2-D"  . 

AND.  HI  .LT. 

0.0)) 

GO 

TO 

50 

K = 

K*  1 

S WRITE  (6  , 

904) 

SO 

IF 

(TWTiN 

F . G t . 0.0) 

GO  TO  60 

K= 

K*  1 

« WRITE  (6, 

905) 

RO 

IF 

(K  .Ew 

. 0)  RETURN 

RETURN  I 

...FORMAT  STATEMENTS 
C 

900  FORMAT (//, 30X," IDENTICAL  INPUT  li NES:  DELX=0.0 — BOUNDARY  ", 
i"L  AY  5ii  COMPUTATIONS  TERMINATED") 

9)1  rO:i(1AT  (//,27X, "NEGATIVE  OR  ZERO  MACH  NUMBER  ENTERED", 
i"--BOUNDA  RY  LAYER  COMPUTATIONS  1LRM1NATED") 

902  FORMAT (//, 19X," NEGATIVE  MSTAR  OR  MSI AR  .GT.  SQ8T ( (G ♦ 1 . ) / (G- 1 . ) ) " 
i "E  NT  ERLD--0GU  ND  Al.  Y LAYER  COMPUTATIONS  TERMINATED") 

9JJ  FORMAT  (//, 1 0X , "NC  (OR  NEGATIVE)  BODY  RADIUS  tNTERED  FOR  AXISYM", 
A"M FTRIC  GE0M2T  KY  — 30UNDARY  LAYER  COMPUTATIONS  TERMINATED”) 

909  FORMAT  (//,  17X,  "MuGATI  VE  CENTERLINE  DISTANCE  ENTERED  FOR  PLANE  ", 
i"2-D  GEOMETRY --BOUNDARY  LAYER  COMPUTATIONS  TERMINATED") 

9 Vi  'ORM  AT  (//,33X,  "NEGATIVE  VALUE  OF  T «T  I NF  ENTERED — BOUNDARY  ", 
.•■"LAYER  COMPUTATIONS  TERMINATED") 

END 


IER 

10 

IER 

20 

IER 

30 

IER 

43 

IER 

50 

IER 

60 

IER 

70 

IER 

80 

IER 

90 

IER 

100 

IER 

113 

IER 

120 

IER 

130 

IER 

140 

IER 

150 

IER 

160 

IER 

170 

IER 

180 

IER 

190 

IER 

200 

IER 

210 

IER 

22  0 

IER 

230 

IER 

240 

IER 

250 

IER 

263 

IER 

273 

IER 

283 

IER 

29  0 

IER 

330 

IER 

310 

IER 

323 

IER 

3 J 3 

IER 

340 

IEH 

353 
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APPENDIX  B 


ERROR  MESSAGES  AND  DISCUSSION 


Program  COMPBL  consists  of  a main  program,  in  which  most 
of  the  calculations  are  performed,  and  five  subprograms:  func- 
tion FJP  evaluates  the  quantity  J'  in  Felsch's  [13]  turbulent 
dissipation  law;  subroutine  TRANS  locates  the  laminar  instabil- 
ity and  transition  points;  subroutine  MSORT  sorts  and  converts 
M and  M::  velocity  input  data;  subroutine  TSORT  sorts  and  con- 
verts 9 and  Tw/Too  temperature  input  data;  and  subroutine  I ERROR 
checks  for  input  errors. 

A listing  of  the  error  messages  generated  by  program  COMPBL 
and  subroutine  I ERROR  is  given  below,  together  with  an  explana- 
tion of  each.  Most  of  these  messages  involve  incorrect  input 
values . 

1.  INPUT  ERROR:  FLOW  MUST  EQUAL  "LAMINAR"  OR  "TURBULENT" 

Literal  variable  FLOW,  entered  in  A10  format,  must  have 
one  of  the  values  given. 

2.  INPUT  ERROR:  GEOM  MUST  EQUAL  "PLANE  2-D"  OR  "ASIXYM" 

Literal  variable  GEOM,  entered  in  A10  format,  must  have 
one  of  the  values  given. 

3.  INPUT  ERROR:  MTYPE  MUST  EQUAL  "MACH"  OR  "MSTAR" 

Literal  variable  MTYPE,  entered  in  A10  format,  must  have 
one  of  the  values  given. 

4.  INPUT  ERROR:  TTYPE  MUST  EQUAL  "THETA"  OR  "TWTINF" 

Literal  variable  TTYPE,  entered  in  A10  format,  must  have 
one  cf  the  values  given. 

5.  INPUT  ERROR:  HSTART  MUST  BE  .GT.  1.515  AND  . LE  . 2.0 

Variable  HS TART , entered  in  NAMELIST  BL,  must  satisfy 
the  inequalities:  1.515  < HSTART  <_  2.0.  (These  inequali- 
ties must  be  satisfied  for  H at  any  location.) 

6.  INPUT  ERROR:  BSTART  MUST  BE.GE.  -0.1988  AND  . LE . 2.0 

Variable  BSTART,  entered  in  NAMELIST  BL,  must  satisfy 
the  following  inequalities:  -0.1988  £ BSTART  <_  2.0. 

7.  INPUT  ERROR:  ZSTART,  FSTINT,  W MUST  BE.GE.  0.0  AND  MINF, 
REINFL,  RL,  RT,  SL,  ST,  EPS  MUST  BE  .GT.  0.0 

Variables  ZSTART,  FSTINT,  W,  MINF,  REINFL,  RL,  RT,  SL, 

ST,  and  EPS,  all  entered  in  NAMELIST  BL,  must  satisfy 
the  following  inequalities:  ZSTART,  FSTINT,  W >_  0.0, 

MINF,  REINFL,  RL,  RT,  SL,  ST,  EPS  > 0.0. 
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8. 


INPUT  ERROR:  G MUST  BE  .GE.  1.0  AND  . LE . 1.67 

Variable  G,  entered  in  NAMELIST  BL,  must  satisfy  the 
following  inequalities:  1.0  £ G £ 1.67. 

9.  INPUT  ERROR:  MINF”  MUST  BE  . LE . SQRT  CCG+1 . )/ (G- 1 . )) 

If  MTYPE  = MSTAR,  variable  MINF,  which  is  entered  in 
NAMELIST  BL,  must  satisfy  the  following  inequality: 

MINF  < /(G+l) 

V CG-1) 

10.  EXECUTION  DELETED 

If  any  of  the  errors  above  are  detected,  this  message 
is  also  printed. 

11.  IDENTICAL  INPUT  LINES:  DELX  = 0 . 0--BOUNDARY  LAYER  COMPU- 
TATIONS TERMINATED 

If  the  X,  Y coordinates  of  successive  points  are  iden- 
tical, the  distance  between  them,  given  by  Eq.  (33)  is 
zero,  which  would  lead  to  division  by  zero. 

12.  NEGATIVE  OR  ZERO  MACH  NUMBER  ENTERED--BOUNDARY  LAYER  COM- 
PUTATIONS TERMINATED 

The  Mach  number  velocity  input  data  must  satisfy  the 
inequality  M >0  (first  point  not  checked  in  order 
to  allow  stagnation  there) . 

13.  NEGATIVE  MSTAR  OR  MSTAR  .GT.  SQRT  ( (G+ 1 . )/ (G- 1 . ) ) ENTERED-- 
BOUNDARY  LAYER  COMPUTATIONS  TERMINATED 

The  M”  velocity  input  data  must  satisfy  the  inequali- 
ites:  0 < M”  _<  /”(G+ 1 )/ CG-  1 ) (first  point  not  checked  in 

order  to  allow  stagnation  there) . 

14.  NO  (OR  NEGATIVE)  BODY  RADIUS  ENTERED  FOR  AXISYMMETRIC 
GEOMETRY --BOUNDARY  LAYER  COMPUTATIONS  TERMINATED 

For  axisymmetric  geometries,  a radius  value,  R,  must 
be  entered,  and  for  all  but  the  first  point  it  must 
satisfy:  R > 0. 

15.  NEGATIVE  CENTERLINE  DISTANCE  ENTERED  FOR  PLANE  2-D  GEOMETRY-- 
BOUNDARY  LAYER  COMPUTATIONS  TERMINATED 

For  plane  geometries  the  centerline-to-boundary  distance, 
R,  must  satisfy  the  inequality:  R ^ 0. 

16.  NEGATIVE  VALUE  OF  TWTINF  ENTERED--BOUNDARY  LAYER  COMPUTA- 
TIONS TERMINATED 

For  TTYPE  = TWTINF,  the  temperature  input  data  must 
satisfy  the  inequality:  TWTINF  >_  0. 

17.  PROBABLE  CONVERGENCE  PROBLEMS:  H EXCEEDS  2 . 0--BOUNDARY 
LAYER  COMPUTATIONS  TERMINATED 
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The  largest  value  that  H can  obtain  is  2.0.  If  H ex- 
ceeds this  value,  it  is  probable  that  the  iterations  for 
H have  become  unstable.  This  may  occur  in  starting 
laminar  boundary  layers  if  an  incorrect  value  of  HSTART 
is  specified  or  if  the  base  points  are  too  widely  sepa- 
rated. Re-examination  of  HSTART,  finer  base  point 
spacing,  and/or  use  of  the  ERROR  = .TRUE,  option  are 
recommended . 

18.  CONVERGENCE  PROBLEMS  FOR  H GIVEN  HSTAR:  EXECUTION  HALTED 

An  iterative  technique  is  utilized  to  solve  for  II  given 
the  corresponding  H”  value.  If  the  required  number  of 
iterations  exceeds  20,  the  message  above  is  written. 

This  error  condition  has  never  been  encountered. 

19.  CONVERGENCE  PROBLEMS  FOR  H IN  MAIN  ITERATION  L00P:EXE- 
CUTION  HALTED 

This  message  is  printed  if  the  number  of  iterations  re- 
quired to  determine  H in  the  iterative  technique  ex- 
plained in  Section  II-A  exceeds  20.  This  error  condi- 
tion has  been  encountered,  although  very  infrequently. 
Increasing  the  number  of  base  points  and/or  use  of 
the  ERROR  = .TRUE,  option  to  determine  the  source  of 
the  problem  are  recommended . 

20.  LAMINAR  SEPARATI 0N--B0UNDARY  LAYER  COMPUTATIONS  TERMINATED 

or 

TURBULENT  SEPARATI 0N--B0UNDARY  LAYER  COMPUTATIONS  TERMINATED 
Although  these  are  not  normally  error  messages,  they  may 
indicate  instability  in  the  iterations  for  H if  they  occur 
at  totally  unexpected  locations,  e.g. , regions  of  favor- 
able pressure  gradient.  This  is  true  because  separa- 
tion is  defined  as:  H £ 1.515.  As  explained  above  in 
Error  17,  the  most  likely  occurrence  of  this  condition 
is  in  starting  laminar  boundary  layer  computations.  Re- 
examination of  the  value  of  HSTART,  finer  base  point 
spacing,  and/or  use  of  the  ERROR  = .TRUE,  option  are  re- 
commended . 
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